#### purpose: analyze terry stop data ####
## Seattle data

#### attach packages ####

rm(list=ls())
suppressPackageStartupMessages(
  
  {
    library(dplyr)
    library(zoo)
    library(estimatr)
    library(tidyverse)
    library(haven)
    library(readxl)
    library(ggthemes)
    library(rdrobust)
    library(gridExtra)
    library(lubridate)
    library(extrafont)
    library(tools)
    library(archivist)
    library(grid)
  }
  
)

#### load and clean data ####

# loading stop data 

load("terry.RData")

terry = terry %>% 
  mutate(blm = ifelse(date >= as.Date("2020-05-29"), 1, 0)) %>% 
  mutate(initiate_off = ifelse(call_type == "-" | call_type == "ONVIEW" | 
                    call_type == "SCHEDULED EVENT (RECURRING)" | 
                      call_type == "ALARM CALL (NOT POLICE ALARM)", 1, 0),
         initiate_civ = 
           ifelse(call_type == "911" | 
                    call_type == "TELEPHONE OTHER, NOT 911" | 
                    call_type == "TEXT MESSAGE", 1, 0)) %>% 
  mutate(hit = ifelse(stop_resolution == "Arrest" | 
                        stop_resolution == "Citation / Infraction" | 
                        stop_resolution == "Offense Report" | 
                        stop_resolution == "Referred for Prosecution", 1, 0)) %>%
  mutate(arrest = ifelse(stop_resolution == "Arrest", 1, 0))

# loading 911 call data
 
load(file = "calldata.RData")

call = 
  call %>% 
  mutate(blm = ifelse(date >= as.Date("2020-05-29"), 1, 0)) %>% 
  filter(date <= as.Date("2021-06-01")) %>% # Filter out everything after call center was transferred to civilian control
  mutate(initiate_off = ifelse(call_type == "-" | call_type == "ONVIEW" | 
                                 call_type == "SCHEDULED EVENT (RECURRING)" | 
                                 call_type == "ALARM CALL (NOT POLICE ALARM)" | 
                                 call_type == "PROACTIVE (OFFICER INITIATED)" | 
                                 call_type == "POLICE (VARDA) ALARM", 1, 0),
         initiate_civ = ifelse(call_type == "911" | 
                    call_type == "TELEPHONE OTHER, NOT 911" | 
                    call_type == "TEXT MESSAGE" | 
                    call_type == "HISTORY CALL (RETRO)" | 
                    call_type == "IN PERSON COMPLAINT", 1, 0))

# load geographic data

load("terry_merged_use.Rdata")
load("crime_merged_use.Rdata")
load("calldata_merged_use.Rdata")

terry_merged_use = terry_merged_use %>% 
  mutate(blm = ifelse(date >= as.Date("2020-05-29"), 1, 0)) %>% 
  mutate(initiate_off = ifelse(call_type == "-" | call_type == "ONVIEW" | 
                                 call_type == "SCHEDULED EVENT (RECURRING)" | 
                                 call_type == "ALARM CALL (NOT POLICE ALARM)", 1, 0),
         initiate_civ = 
           ifelse(call_type == "911" | 
                    call_type == "TELEPHONE OTHER, NOT 911" | 
                    call_type == "TEXT MESSAGE", 1, 0))

calldata_merged_use = 
  calldata_merged_use %>% 
  mutate(blm = ifelse(date >= as.Date("2020-05-29"), 1, 0)) %>% 
  filter(date <= as.Date("2021-06-01")) %>% # Filter out everything after call center was transferred to civilian control
  mutate(initiate_off = ifelse(call_type == "-" | call_type == "ONVIEW" | 
                                 call_type == "SCHEDULED EVENT (RECURRING)" | 
                                 call_type == "ALARM CALL (NOT POLICE ALARM)" | 
                                 call_type == "PROACTIVE (OFFICER INITIATED)" | 
                                 call_type == "POLICE (VARDA) ALARM", 1, 0),
         initiate_civ = ifelse(call_type == "911" | 
                                 call_type == "TELEPHONE OTHER, NOT 911" | 
                                 call_type == "TEXT MESSAGE" | 
                                 call_type == "HISTORY CALL (RETRO)" | 
                                 call_type == "IN PERSON COMPLAINT", 1, 0)) %>%
  filter(initiate_off == 1)


# gen time data

call$time_arrived2 = parse_date_time(call$time_arrived, "%m/%d/%Y %H:%M:%S %p")
call$time_queued2 = parse_date_time(call$time_queued, "%m/%d/%Y %H:%M:%S %p")

call$time_arrived2 = 
  as.POSIXct(call$time_arrived, format="%m/%d/%Y %H:%M:%S %p", tz="UTC")
call$time_queued2 = 
  as.POSIXct(call$time_queued, format="%m/%d/%Y %H:%M:%S %p", tz="UTC")
call$response_time = (call$time_arrived2 - call$time_queued2) / 60
the_index = call$response_time < 0
call = call %>% filter(!the_index)

calldata_merged_use$time_arrived2 = parse_date_time(calldata_merged_use$time_arrived, "%m/%d/%Y %H:%M:%S %p")
calldata_merged_use$time_queued2 = parse_date_time(calldata_merged_use$time_queued, "%m/%d/%Y %H:%M:%S %p")

calldata_merged_use$time_arrived2 = 
  as.POSIXct(calldata_merged_use$time_arrived, format="%m/%d/%Y %H:%M:%S %p", tz="UTC")
calldata_merged_use$time_queued2 = 
  as.POSIXct(calldata_merged_use$time_queued, format="%m/%d/%Y %H:%M:%S %p", tz="UTC")
calldata_merged_use$response_time = (calldata_merged_use$time_arrived2 - calldata_merged_use$time_queued2) / 60
the_index = calldata_merged_use$response_time < 0
calldata_merged_use = calldata_merged_use %>% filter(!the_index)

# loading in crime data 

load(file = "crimedata.RData")

crime = 
  crime %>% 
  mutate(date = as.Date(paste(substring(report_date, 7, 10),
                              substring(report_date, 1, 2),
                              substring(report_date, 4, 5), sep = "-"))) %>% 
  mutate(blm = ifelse(date >= as.Date("2020-05-29"), 1, 0)) %>% 
  mutate(cat_society = ifelse(crime_against_category == "SOCIETY", 1, 0),
         cat_property = ifelse(crime_against_category == "PROPERTY", 1, 0),
         cat_person = ifelse(crime_against_category == "PERSON", 1, 0),
         cat_notcrime = ifelse(crime_against_category == "NOT_A_CRIME", 1, 0)) %>% 
  filter(offense != "Identity Theft")


#### analyzing data ####

#### ot plot, blm protest, depolicing (stops + 911) ####

# rdit analysis

terry_ts = terry %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) 

call_ts = call %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            response_time = mean(response_time, na.rm = TRUE)) 

call_ts2 = call %>% 
  mutate(count = 1) %>% 
  filter(initiate_civ == 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            response_time = mean(response_time, na.rm = TRUE))

call_ts3 = call %>% 
  mutate(count = 1) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            response_time = mean(response_time, na.rm = TRUE))

terry_ts$blmrv = 
  terry_ts$date - min(terry_ts$date[terry_ts$blm == 1])
call_ts$blmrv = 
  call_ts$date - min(call_ts$date[call_ts$blm == 1])
call_ts2$blmrv = 
  call_ts2$date - min(call_ts2$date[call_ts2$blm == 1])
call_ts3$blmrv = 
  call_ts3$date - min(call_ts3$date[call_ts3$blm == 1])

rdout = 
  with(terry_ts, rdrobust(x = blmrv, y = count, p = 1, kernel = "tri"))
rdout2 = 
  with(call_ts, rdrobust(x = blmrv, y = count, p = 1, kernel = "tri"))
rdout3 = 
  with(call_ts2, rdrobust(x = blmrv, y = response_time, p = 1, kernel = "tri"))

# plot 

depolp1 = terry %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>% 
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .5) + 
  labs(x = "Date", y = "Count",
       title = "A. Terry Stops") + 
  annotate("text", label = paste0("RDiT Coef: ",
                                  round(rdout$coef[1]),
                                  "\nSE: ",
                                  round(rdout$se[3]),
                                  "\np < 0.001",
                                  "\npre-BLM SD: ",
                                  round(sd(terry_ts$count[terry_ts$blm == 0]))),
           x = as.Date("2021-03-01"),
           y = 40,
           family = "serif",
           size = 1.5) + 
  theme_tufte(base_size = 9)

depolp2 = call %>% 
  mutate(count = 1) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>% 
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .5) + 
  labs(x = "Date", y = "Count",
       title = "B. Officer-Initiated Calls") + 
  annotate("text", label = paste0("RDiT Coef: ",
                                  round(rdout2$coef[1]),
                                  "\nSE: ",
                                  round(rdout2$se[3]),
                                  "\np < 0.001",
                                  "\npre-BLM SD: ",
                                  round(sd(call_ts$count[call_ts$blm == 0]))),
           x = as.Date("2021-03-01"),
           y = 500,
           family = "serif",
           size = 1.5) + 
  theme_tufte(base_size = 9)


depolp3 = call %>% 
  mutate(count = 1) %>% 
  filter(initiate_civ == 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            response_time = mean(response_time, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>% 
  ggplot() + 
  geom_point(aes(x = date, y = response_time),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = response_time, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             ize = .5) + 
  annotate("text", label = paste0("RDiT Coef: ",
                                  round(rdout3$coef[1]),
                                  "\nSE: ",
                                  round(rdout3$se[3]),
                                  "\np < 0.001",
                                  "\npre-BLM SD: ",
                                  round(sd(call_ts2$response_time[call_ts$blm == 0]))),
           x = as.Date("2021-03-01"),
           y = 175,
           family = "serif",
           size = 1.5) + 
  labs(x = "Date", y = "Mean Response Time",
       title = "C. Response Times") + 
  ylim(0, 200) + 
  theme_tufte(base_size = 9)


depolpGrob = arrangeGrob(depolp1, depolp2, depolp3, ncol = 3)
ggsave(plot = depolpGrob, filename = "depol.png", width = 8, height = 2)

#### ot plot, blm protest, terry stops depolicing disagg ####

disaggterry1 = terry %>% 
  mutate(count = 1) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>% 
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .2,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .5,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2) + 
  labs(x = "Date", y = "Count",
       title = "A. Terry Stops (Officer-Initiated)") + 
  theme_tufte()

disaggterry2 = terry %>% 
  mutate(count = 1) %>% 
  filter(initiate_civ == 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>% 
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .2,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .5,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2) + 
  labs(x = "Date", y = "Count",
       title = "B. Terry Stops (Civilian-Initiated)") + 
  theme_tufte()

disaggterryGrob = arrangeGrob(disaggterry1, disaggterry2, ncol = 2)
ggsave(plot = disaggterryGrob, filename = "terrydisagg.png", width = 8, height = 2.5)

#### ot plot, effectiveness #### 

# rdit arrest

terry_arr_ts = 
  terry %>% 
  mutate(count = 1,
         weapnrep = ifelse(weapon_type == "None" | weapon_type == "-", 1, 0)) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            weapnrep = mean(weapnrep, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            arrest = mean(arrest, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2019-06-01"))
prop.table(table(terry$weapon_type))

mean((1 - terry_arr_ts$weapnrep)[terry_arr_ts$blm == 0])
terry_arr_ts %>% 
  ggplot() + 
  geom_point(aes(x = date, y = weapnrep)) + 
  geom_smooth(aes(x = date, y = weapnrep, group = blm))


terry_arr_ts$blmrv = 
  terry_arr_ts$date - min(terry_arr_ts$date[terry_arr_ts$blm == 1])

rdout = 
  with(terry_arr_ts, rdrobust(x = blmrv, y = arrest, p = 1, kernel = "tri", all = TRUE))

# rdit hit rate

terry_hit_ts = terry %>% 
  mutate(count = 1) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  summarize(hit = mean(hit, na.rm = TRUE),
            count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>%
  filter(date >= as.Date("2019-06-01"))


terry_hit_ts$blmrv = 
  terry_hit_ts$date - min(terry_hit_ts$date[terry_hit_ts$blm == 1])

rdout2 = 
  with(terry_hit_ts, rdrobust(x = blmrv, y = hit, p = 1, kernel = "tri", all = TRUE))

arrest_rate = 
  terry %>% 
  mutate(count = 1) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            arrest = mean(arrest, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2019-06-01")) %>% 
  ggplot() + 
  geom_point(aes(x = date, y = arrest),
             size = .2,
             alpha = .2) +
  geom_smooth(aes(x = date, y = arrest, group = blm),
              se = FALSE,
              size = .5,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2) + 
  labs(x = "Date", y = "Rate",
       title = "A. Arrest Rate") + 
  annotate("text", label = paste0("RDiT Coef: ",
                                  round(rdout$coef[1], 2),
                                  "\nSE: ",
                                  round(rdout$se[3], 2),
                                  "\np < 0.05",
                                  "\npre-BLM SD: ",
                                  round(sd(terry_arr_ts$arrest[terry_arr_ts$blm == 0]), 2)),
           x = as.Date("2021-03-01"),
           y = .75,
           family = "serif",
           size = 2.25) + 
  theme_tufte()

hit_rate = terry_hit_ts %>% 
  ggplot() + 
  geom_point(aes(x = date, y = hit),
             size = .2,
             alpha = .2) +
  geom_smooth(aes(x = date, y = hit, group = blm),
              se = FALSE,
              size = .5,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2) + 
  labs(x = "Date", y = "Rate",
       title = "B. Hit Rate") + 
  annotate("text", label = paste0("RDiT Coef: ",
                                  round(rdout2$coef[1], 2),
                                  "\nSE: ",
                                  round(rdout2$se[3], 2),
                                  "\np < 0.05",
                                  "\npre-BLM SD: ",
                                  round(sd(terry_hit_ts$hit[terry_hit_ts$blm == 0]), 2)),
           x = as.Date("2021-03-01"),
           y = .75,
           family = "serif",
           size = 2.25) + 
  theme_tufte()

grobEff = arrangeGrob(arrest_rate, hit_rate, ncol = 2)
ggsave(plot = grobEff, filename = "eff.png", width = 8, height = 2.5)

#### ot plot, effectiveness, assessing raw #### 

# arrests

terry_arr_ts = 
  terry %>% 
  mutate(count = 1) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            arrest = sum(arrest, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2019-06-01"))


terry_arr_ts$blmrv = 
  terry_arr_ts$date - min(terry_arr_ts$date[terry_arr_ts$blm == 1])

rdout = 
  with(terry_arr_ts, rdrobust(x = blmrv, y = arrest, p = 1, kernel = "tri", all = TRUE))

# hits

terry_hit_ts = terry %>% 
  mutate(count = 1) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  summarize(hit = sum(hit, na.rm = TRUE),
            count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>%
  filter(date >= as.Date("2019-06-01"))

terry_hit_ts$blmrv = 
  terry_hit_ts$date - min(terry_hit_ts$date[terry_hit_ts$blm == 1])

rdout2 = 
  with(terry_hit_ts, rdrobust(x = blmrv, y = hit, p = 1, kernel = "tri", all = TRUE))

arrest_num = 
  terry %>% 
  mutate(count = 1) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            arrest = sum(arrest, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2019-06-01")) %>% 
  ggplot() + 
  geom_point(aes(x = date, y = arrest),
             size = .2,
             alpha = .2) +
  geom_smooth(aes(x = date, y = arrest, group = blm),
              se = FALSE,
              size = .5,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2) + 
  labs(x = "Date", y = "Count",
       title = "A. Arrests") + 
  annotate("text", label = paste0("RDiT Coef: ",
                                  round(rdout$coef[1], 2),
                                  "\nSE: ",
                                  round(rdout$se[3], 2),
                                  "\np < 0.05",
                                  "\npre-BLM SD: ",
                                  round(sd(terry_arr_ts$arrest[terry_arr_ts$blm == 0]), 2)),
           x = as.Date("2021-03-01"),
           y = 6,
           family = "serif",
           size = 2.25) + 
  theme_tufte()

hit_num = terry_hit_ts %>% 
  ggplot() + 
  geom_point(aes(x = date, y = hit),
             size = .2,
             alpha = .2) +
  geom_smooth(aes(x = date, y = hit, group = blm),
              se = FALSE,
              size = .5,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2) + 
  labs(x = "Date", y = "Count",
       title = "B. Hits") + 
  annotate("text", label = paste0("RDiT Coef: ",
                                  round(rdout2$coef[1], 2),
                                  "\nSE: ",
                                  round(rdout2$se[3], 2),
                                  "\np < 0.001",
                                  "\npre-BLM SD: ",
                                  round(sd(terry_hit_ts$hit[terry_hit_ts$blm == 0]), 2)),
           x = as.Date("2021-03-01"),
           y = 6,
           family = "serif",
           size = 2.25) + 
  theme_tufte()

grobEff = arrangeGrob(arrest_num, hit_num, ncol = 2)
ggsave(plot = grobEff, filename = "eff2.png", width = 8, height = 2.5)

#### ot plot, disparities #### 

terry_arr_ts = 
  terry %>% 
  mutate(count = 1,
         r_unk = ifelse(subject_race == "Unknown" | subject_race == "-", 1, 0),
         r_wht = ifelse(subject_race == "White", 1, 0),
         r_asn = ifelse(subject_race == "Asian", 1, 0),
         r_blk = ifelse(subject_race == "Black or African American", 1, 0)) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            arrest = mean(arrest, na.rm = TRUE),
            r_wht = sum(r_wht, na.rm = TRUE),
            r_blk = sum(r_blk, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2019-06-01")) %>% 
  mutate(wht_rate = (r_wht / (.673 * 737015)) * 10000,
         blk_rate = (r_blk / (.073 * 737015)) * 10000) %>% 
  mutate(rr_bw = blk_rate / wht_rate) %>% 
  mutate(rr_bw = ifelse(rr_bw == Inf | is.na(rr_bw), NA, rr_bw))


# running variable

terry_arr_ts$blmrv = 
  terry_arr_ts$date - min(terry_arr_ts$date[terry_arr_ts$blm == 1])

# rdit 

rdout = 
  with(terry_arr_ts, rdrobust(x = blmrv, y = rr_bw, p = 1, kernel = "tri", all = TRUE))

# now, plot 

stoprp2 = terry_arr_ts %>% 
  ggplot() + 
  geom_point(aes(x = date, y = rr_bw),
             size = .2,
             alpha = .2) +
  geom_smooth(aes(x = date, y = rr_bw, group = blm),
              se = FALSE,
              size = .5,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2) + 
  labs(x = "Date", y = "Black/White Stop Rate Ratio",
       title = "B. Black/White Stop Rate Ratio") + 
  annotate("text", label = paste0("RDiT Coef: ",
                                  round(rdout$coef[1], 2),
                                  "\nSE: ",
                                  round(rdout$se[3], 2),
                                  "\np < 0.01",
                                  "\npre-BLM SD: ",
                                  round(sd(terry_arr_ts$rr_bw[terry_arr_ts$blm == 0], na.rm = TRUE), 2)),
           x = as.Date("2021-03-01"),
           y = 40,
           family = "serif",
           size = 2.25) + 
  theme_tufte(base_size = 9)

stoprp1 = terry_arr_ts %>% 
  ggplot() + 
  geom_point(aes(x = date, y = blk_rate),
             size = .2,
             alpha = .2) +
  geom_smooth(aes(x = date, y = blk_rate, group = blm),
              se = FALSE,
              size = .5,
              col = "black") + 
  geom_point(aes(x = date, y = wht_rate),
             size = .2,
             alpha = .2,
             col = "grey") +
  geom_smooth(aes(x = date, y = wht_rate, group = blm),
              se = FALSE,
              size = .5,
              col = "grey") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2) + 
  labs(x = "Date", y = "Stop Rate (per 10,000)",
       title = "A. Black (Black) and White (Grey) Stop Rate") + 
  theme_tufte(base_size = 9)

stoprp_grob = arrangeGrob(stoprp1, stoprp2, ncol = 2)
ggsave(filename = "stoprates.png", width = 8, height = 2, plot = stoprp_grob)

#### ot plot, crime, assessing raw #### 

# time series dataset

crime_ts = 
  crime %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(cat_property = sum(cat_property, na.rm = TRUE),
            cat_person = sum(cat_person, na.rm = TRUE),
            cat_society = sum(cat_society, na.rm = TRUE),
            count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>%
  filter(date >= as.Date("2016-01-01"))

crime_ts$blmrv = crime_ts$date - min(crime_ts$date[crime_ts$blm == 1])

# rdits 

rdit_crime1 = with(crime_ts, rdrobust(x = blmrv, y = count, p = 1, all = TRUE))
rdit_crime2 = with(crime_ts, rdrobust(x = blmrv, y = cat_property, p = 1, all = TRUE))
rdit_crime3 = with(crime_ts, rdrobust(x = blmrv, y = cat_society, p = 1, all = TRUE))
rdit_crime4 = with(crime_ts, rdrobust(x = blmrv, y = cat_person, p = 1, all = TRUE))

# plot 


crimep1 = crime_ts %>% 
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .2,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .5,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2) + 
  labs(x = "Date", y = "Count",
       title = "A. All") + 
  annotate("text", label = paste0("RDiT Coef: ",
                                  round(rdit_crime1$coef[1], 2),
                                  "\nSE: ",
                                  round(rdit_crime1$se[3], 2),
                                  "\np < 0.001",
                                  "\npre-BLM SD: ",
                                  round(sd(crime_ts$count[crime_ts$blm == 0]), 2)),
           x = as.Date("2021-03-01"),
           y = 600,
           family = "serif",
           size = 1) + 
  theme_tufte(base_size = 8)

crimep2 = crime_ts %>% 
  ggplot() + 
  geom_point(aes(x = date, y = cat_property),
             size = .2,
             alpha = .2) +
  geom_smooth(aes(x = date, y = cat_property, group = blm),
              se = FALSE,
              size = .5,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2) + 
  labs(x = "Date", y = "Count",
       title = "B. Property") + 
  annotate("text", label = paste0("RDiT Coef: ",
                                  round(rdit_crime2$coef[1], 2),
                                  "\nSE: ",
                                  round(rdit_crime2$se[3], 2),
                                  "\np < 0.001",
                                  "\npre-BLM SD: ",
                                  round(sd(crime_ts$cat_property[crime_ts$blm == 0]), 2)),
           x = as.Date("2021-03-01"),
           y = 500,
           family = "serif",
           size = 1) + 
  theme_tufte(base_size = 8)

crimep3 = crime_ts %>% 
  ggplot() + 
  geom_point(aes(x = date, y = cat_society),
             size = .2,
             alpha = .2) +
  geom_smooth(aes(x = date, y = cat_society, group = blm),
              se = FALSE,
              size = .5,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2) + 
  labs(x = "Date", y = "Count",
       title = "C. Against Society") + 
  annotate("text", label = paste0("RDiT Coef: ",
                                  round(rdit_crime3$coef[1], 2),
                                  "\nSE: ",
                                  round(rdit_crime3$se[3], 2),
                                  "\np < 0.001",
                                  "\npre-BLM SD: ",
                                  round(sd(crime_ts$cat_society[crime_ts$blm == 0]), 2)),
           x = as.Date("2021-03-01"),
           y = 40,
           family = "serif",
           size = 1) + 
  theme_tufte(base_size = 8)

crimep4 = crime_ts %>% 
  ggplot() + 
  geom_point(aes(x = date, y = cat_person),
             size = .2,
             alpha = .2) +
  geom_smooth(aes(x = date, y = cat_person, group = blm),
              se = FALSE,
              size = .5,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2) + 
  labs(x = "Date", y = "Count",
       title = "D. Against Person") + 
  annotate("text", label = paste0("RDiT Coef: ",
                                  round(rdit_crime4$coef[1], 2),
                                  "\nSE: ",
                                  round(rdit_crime4$se[3], 2),
                                  "\np < ", round(rdit_crime4$pv[3], 2),
                                  "\npre-BLM SD: ",
                                  round(sd(crime_ts$cat_person[crime_ts$blm == 0]), 2)),
           x = as.Date("2021-03-01"),
           y = 55,
           family = "serif",
           size = 1) + 
  theme_tufte(base_size = 8)

crime_grob = arrangeGrob(crimep1, crimep2, crimep3, crimep4, ncol = 4)
ggsave(filename = "crimedat.png", plot = crime_grob, width = 8, height = 2)

#### analyses --- depolicing #### 

# datasets

datasets_depol = list(terry_ts %>% as.data.frame %>% mutate(count_terry = count) %>% filter(date >= as.Date("2016-01-01")),
                      call_ts %>% as.data.frame %>% mutate(count_calls = count) %>% filter(date >= as.Date("2016-01-01")),
                      call_ts %>% as.data.frame %>% filter(date >= as.Date("2016-01-01")))

# outcomes

outcomes = c("count_terry", "count_calls", "response_time")

# polynomials

polys = c(0, 1, 2)

# kernels 

kerns = c("tri", "uni", "epa")

# loop for analysis 

dataset_list = as.list(rep(NA, length(datasets_depol)))

for (i in 1:length(datasets_depol)) {
  
    poly_list = as.list(rep(NA, length(polys)))
    
    for (j in 1:length(polys)) {
      
      print(paste0("J Iteration ", j))
      
      kern_list = matrix(nrow = length(kerns), ncol = 6)
      
      for (l in 1:length(kerns)) {
        
        print(paste0("L Iteration ", l))
        
        outrd = 
          rdrobust(x = as.numeric(datasets_depol[[i]][, "blmrv"]),
                                  y = datasets_depol[[i]][, outcomes[i]],
                                  p = polys[j],
                                  kernel = kerns[l])
        
        kern_list[l, 1] = outrd$coef[3]
        kern_list[l, 2] = outrd$se[3]
        kern_list[l, 3] = outrd$pv[3]
        kern_list[l, 4] = polys[j]
        kern_list[l, 5] = kerns[l]
        kern_list[l, 6] = outcomes[i]
        
      }
      
      poly_list[[j]] = kern_list %>% 
        as.data.frame %>% 
        `colnames<-` (c("est", "se", "pv", "polys", "kerns", "outcomes")) %>% 
        mutate(est = as.numeric(as.character(est)),
               se = as.numeric(as.character(se)),
               pv = as.numeric(as.character(pv)),
               polys = as.numeric(as.character(polys)),
               kerns = as.character(kerns),
               outcomes = as.character(outcomes))
      
    }
  
  dataset_list[[i]] = poly_list
  
}

depol_rdit_res = 
  dataset_list %>% 
  unlist(recursive = FALSE) %>% 
  do.call(rbind.data.frame, .) %>% 
  mutate(outcomes = factor(outcomes,
                           levels = c("count_terry",
                                     "count_calls",
                                     "response_time"),
                           labels = c("A. Terry Stops",
                                      "B. 911 Calls",
                                      "C. Response Time")),
         polys = factor(polys, levels = c(0, 1, 2),
                        labels = c("DIM", "Linear", "Quadratic")),
         kerns = factor(kerns, levels = c("tri", "uni", "epa"),
                        labels = c("Triangular",
                                   "Uniform",
                                   "Epanechnikov")))
  
depol_rdit_res %>% 
  ggplot() + 
  geom_point(aes(x = polys, y = est, col = kerns),
             position = position_dodge(.4)) +
  geom_errorbar(aes(x = polys, 
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    col = kerns),
                position = position_dodge(.4),
                width = 0,
                size = .2) + 
  geom_errorbar(aes(x = polys, 
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    col = kerns),
                position = position_dodge(.4),
                width = 0,
                size = .4) + 
  geom_hline(yintercept = 0) + 
  scale_color_grey(start = 0, end = .6) + 
  facet_wrap(~ outcomes, scales = "free") + 
  labs(x = "Running Variable Specification",
       y = "Coefficient\n(BLM Protest)",
       col = "Kernel") + 
  theme_tufte()

ggsave(plot = last_plot(), width = 8, height = 2, filename = "resrdit1.png")

#### analyses --- depolicing (closer bandwidth) #### 

# datasets

datasets_depol = list(terry_ts %>% as.data.frame %>% mutate(count_terry = count) %>% filter(date >= as.Date("2016-01-01")),
                      call_ts %>% as.data.frame %>% mutate(count_calls = count) %>% filter(date >= as.Date("2016-01-01")),
                      call_ts %>% as.data.frame %>% filter(date >= as.Date("2016-01-01")),
                      call_ts2 %>% as.data.frame %>% mutate(count_calls_civ = count) %>% filter(date >= as.Date("2016-01-01")),
                      call_ts3 %>% as.data.frame %>% mutate(count_calls_off = count) %>% filter(date >= as.Date("2016-01-01")))

# outcomes

outcomes = c("count_terry", "count_calls", "response_time", "count_calls_civ", "count_calls_off")

# polynomials

polys = c(0, 1, 2)

# kernels 

kerns = c("tri", "uni", "epa")

# loop for analysis 

dataset_list = as.list(rep(NA, length(datasets_depol)))

for (i in 1:length(datasets_depol)) {
  
  poly_list = as.list(rep(NA, length(polys)))
  
  for (j in 1:length(polys)) {
    
    print(paste0("J Iteration ", j))
    
    kern_list = as.list(rep(NA, length(kerns)))
    
    for (l in 1:length(kerns)) {
      
      print(paste0("L Iteration ", l))
     
      mat_bw = matrix(NA, 100, 7)
      
      for(z in 10:100) {
        
        print(paste0("Z Iteration ", z))
        
        outrd = 
          rdrobust(x = as.numeric(datasets_depol[[i]][, "blmrv"]),
                   y = datasets_depol[[i]][, outcomes[i]],
                   p = polys[j],
                   h = z,
                   kernel = kerns[l])
        
        mat_bw[z, 1] = outrd$coef[3]
        mat_bw[z, 2] = outrd$se[3]
        mat_bw[z, 3] = outrd$pv[3]
        mat_bw[z, 4] = polys[j]
        mat_bw[z, 5] = kerns[l]
        mat_bw[z, 6] = outcomes[i]
        mat_bw[z, 7] = z

      }
      
      kern_list[[l]] = mat_bw %>%
        as.data.frame() %>%
        `colnames<-` (c("Coefficient", "SE", "Pvalue", "Polynomial", "Kernel", "Outcome", "BW"))%>% 
        mutate(Coefficient = as.numeric(as.character(Coefficient)),
               SE = as.numeric(as.character(SE)),
               Pvalue = as.numeric(as.character(Pvalue)),
               Polynomial = as.numeric(as.character(Polynomial)),
               Kernel = as.character(Kernel),
               Outcome = as.character(Outcome),
               BW = as.character(BW))
    }
    
    poly_list[[j]] = kern_list
    
  }
  
  dataset_list[[i]] = poly_list
  
}

depol_rdit_res = dataset_list %>% 
    unlist(recursive = FALSE) %>%
      unlist(recursive = FALSE) %>%
      do.call(rbind.data.frame, .) %>% 
          mutate(Outcome = factor(Outcome,
                           levels = c("count_terry",
                                      "count_calls",
                                      "response_time",
                                      "count_calls_civ",
                                      "count_calls_off"),
                           labels = c("A. Terry Stops",
                                      "B. 911 Calls Total",
                                      "C. Response Time",
                                      "D. Civilian Initiated Calls",
                                      "E. Officer Initiated Calls")),
              Polynomial = factor(Polynomial, levels = c(0, 1, 2),
                        labels = c("DIM", "Linear", "Quadratic")),
              Kernel = factor(Kernel, levels = c("tri", "uni", "epa"),
                        labels = c("Triangular",
                                   "Uniform",
                                   "Epanechnikov"))) %>% 
            filter(!is.na(BW))

depol_rdit_res = depol_rdit_res %>% 
  mutate(kern_poly = factor(paste0(Kernel, ", ", Polynomial), levels = c("Triangular, DIM", 
                                                                         "Triangular, Linear",
                                                                         "Triangular, Quadratic",
                                                                         "Uniform, DIM",
                                                                         "Uniform, Linear",
                                                                         "Uniform, Quadratic",
                                                                          "Epanechnikov, DIM",
                                                                          "Epanechnikov, Linear",
                                                                          "Epanechnikov, Quadratic")))

depol_rdit_res$BW = as.numeric(depol_rdit_res$BW)

depol_rdit_res %>% filter(Outcome == "A. Terry Stops") %>% 
      ggplot + geom_point(aes(x = BW, y = Coefficient)) + 
        geom_errorbar(aes(x = BW, ymin = Coefficient - 1.96*SE, ymax = Coefficient + 1.96*SE), width = 0) + 
        facet_wrap(~kern_poly, ncol = 3) + 
        geom_hline(yintercept = 0) + 
        theme_tufte() +
        labs(x = "Bandwidth", y = "RDIT Coefficient (BLM)", title = "Terry stops") +
        scale_x_continuous(breaks = seq(from = 10, to = 100, by = 10))

ggsave(plot = last_plot(), filename = "bw_plot_stops.png", width = 8, height = 6)

depol_rdit_res %>% filter(Outcome == "B. 911 Calls Total") %>% 
  ggplot + geom_point(aes(x = BW, y = Coefficient)) + 
  geom_errorbar(aes(x = BW, ymin = Coefficient - 1.96*SE, ymax = Coefficient + 1.96*SE), width = 0) + 
  facet_wrap(~kern_poly, ncol = 3) + 
  geom_hline(yintercept = 0) + 
  theme_tufte() +
  labs(x = "Bandwidth", y = "RDIT Coefficient (BLM)", title = "911 calls (total)") +
  scale_x_continuous(breaks = seq(from = 10, to = 100, by = 10))

ggsave(plot = last_plot(), filename = "bw_plot_calls_total.png", width = 8, height = 6)

depol_rdit_res %>% filter(Outcome == "C. Response Time") %>% 
  ggplot + geom_point(aes(x = BW, y = Coefficient)) + 
  geom_errorbar(aes(x = BW, ymin = Coefficient - 1.96*SE, ymax = Coefficient + 1.96*SE), width = 0) + 
  facet_wrap(~kern_poly, ncol = 3) + 
  geom_hline(yintercept = 0) + 
  theme_tufte() +
  labs(x = "Bandwidth", y = "RDIT Coefficient (BLM)", title = "Response time") +
  scale_x_continuous(breaks = seq(from = 10, to = 100, by = 10))

ggsave(plot = last_plot(), filename = "bw_plot_resptime.png", width = 8, height = 6)

depol_rdit_res %>% filter(Outcome == "D. Civilian Initiated Calls") %>% 
  ggplot + geom_point(aes(x = BW, y = Coefficient)) + 
  geom_errorbar(aes(x = BW, ymin = Coefficient - 1.96*SE, ymax = Coefficient + 1.96*SE), width = 0) + 
  facet_wrap(~kern_poly, ncol = 3) + 
  geom_hline(yintercept = 0) + 
  theme_tufte() +
  labs(x = "Bandwidth", y = "RDIT Coefficient (BLM)", title = "Civilian initiated 911 calls") +
  scale_x_continuous(breaks = seq(from = 10, to = 100, by = 10))

ggsave(plot = last_plot(), filename = "bw_plot_calls_civ.png", width = 8, height = 6)

depol_rdit_res %>% filter(Outcome == "E. Officer Initiated Calls") %>% 
  ggplot + geom_point(aes(x = BW, y = Coefficient)) + 
  geom_errorbar(aes(x = BW, ymin = Coefficient - 1.96*SE, ymax = Coefficient + 1.96*SE), width = 0) + 
  facet_wrap(~kern_poly, ncol = 3) + 
  geom_hline(yintercept = 0) + 
  theme_tufte() +
  labs(x = "Bandwidth", y = "RDIT Coefficient (BLM)", title = "Officer initiated 911 calls") +
  scale_x_continuous(breaks = seq(from = 10, to = 100, by = 10))

ggsave(plot = last_plot(), filename = "bw_plot_calls_off.png", width = 8, height = 6)

#### Depolicing -- Plots for Yale paper ####

depol_plot1 = depol_rdit_res %>% filter(Outcome == "E. Officer Initiated Calls") %>%
  filter(Polynomial == "Linear") %>%
  filter(Kernel == "Uniform") %>%
  ggplot + geom_point(aes(x = BW, y = Coefficient)) + 
  geom_errorbar(aes(x = BW, ymin = Coefficient - 1.96*SE, ymax = Coefficient + 1.96*SE), width = 0) + 
  facet_wrap(~kern_poly, ncol = 3) + 
  geom_hline(yintercept = 0) + 
  theme_tufte() +
  labs(x = "Bandwidth", y = "RDIT Coefficient (BLM)", title = "A. Officer initiated 911 calls") +
  scale_x_continuous(breaks = seq(from = 10, to = 100, by = 10))

ggsave(plot = last_plot(), filename = "bw_plot_calls_off_paper.png", width = 8, height = 6)

depol_plot2 = depol_rdit_res %>% filter(Outcome == "A. Terry Stops") %>%
  filter(Polynomial == "Linear") %>%
  filter(Kernel == "Uniform") %>%
  ggplot + geom_point(aes(x = BW, y = Coefficient)) + 
  geom_errorbar(aes(x = BW, ymin = Coefficient - 1.96*SE, ymax = Coefficient + 1.96*SE), width = 0) + 
  facet_wrap(~kern_poly, ncol = 3) + 
  geom_hline(yintercept = 0) + 
  theme_tufte() +
  labs(x = "Bandwidth", y = "RDIT Coefficient (BLM)", title = "B. Terry Stops") +
  scale_x_continuous(breaks = seq(from = 10, to = 100, by = 10))

ggsave(plot = last_plot(), filename = "bw_plot_terry_paper.png", width = 8, height = 6)

depol_plot3 = depol_rdit_res %>% filter(Outcome == "C. Response Time") %>%
  filter(Polynomial == "Linear") %>%
  filter(Kernel == "Uniform") %>%
  ggplot + geom_point(aes(x = BW, y = Coefficient)) + 
  geom_errorbar(aes(x = BW, ymin = Coefficient - 1.96*SE, ymax = Coefficient + 1.96*SE), width = 0) + 
  facet_wrap(~kern_poly, ncol = 3) + 
  geom_hline(yintercept = 0) + 
  theme_tufte() +
  labs(x = "Bandwidth", y = "RDIT Coefficient (BLM)", title = "D. Response Time") +
  scale_x_continuous(breaks = seq(from = 10, to = 100, by = 10))

ggsave(plot = last_plot(), filename = "bw_plot_resp_paper.png", width = 8, height = 6)

depol_grob_paper = arrangeGrob(depol_plot1, depol_plot2, ncol = 2)
ggsave(filename = "depol_uni_lin.png", plot = depol_grob_paper, width = 8, height = 2)

#### analyses --- efficiency #### 

# datasets

terry_arr_ts = 
  terry %>% 
  mutate(count = 1,
         r_unk = ifelse(subject_race == "Unknown" | subject_race == "-", 1, 0)) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            arrest = mean(arrest, na.rm = TRUE),
            r_unk = mean(r_unk, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2019-06-01"))

terry_arr_ts$blmrv = 
  terry_arr_ts$date - min(terry_arr_ts$date[terry_arr_ts$blm == 1])

# rdit hit rate

terry_hit_ts = terry %>% 
  mutate(count = 1) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  mutate(hit = ifelse(stop_resolution == "Arrest" | 
                        stop_resolution == "Citation / Infraction" | 
                        stop_resolution == "Offense Report" | 
                        stop_resolution == "Referred for Prosecution", 1, 0)) %>% 
  summarize(hit = mean(hit, na.rm = TRUE),
            count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>%
  filter(date >= as.Date("2019-06-01"))

terry_hit_ts$blmrv = 
  terry_hit_ts$date - min(terry_hit_ts$date[terry_hit_ts$blm == 1])

datasets_eff = list(terry_arr_ts %>% as.data.frame,
                    terry_hit_ts %>% as.data.frame)

# outcomes

outcomes = c("arrest", "hit")

# polynomials

polys = c(0, 1, 2)

# kernels 

kerns = c("tri", "uni", "epa")

# loop for analysis 

dataset_list = as.list(rep(NA, length(datasets_eff)))

for (i in 1:length(datasets_eff)) {
  
  poly_list = as.list(rep(NA, length(polys)))
  
  for (j in 1:length(polys)) {
    
    print(paste0("J Iteration ", j))
    
    kern_list = matrix(nrow = length(kerns), ncol = 6)
    
    for (l in 1:length(kerns)) {
      
      print(paste0("L Iteration ", l))
      
      outrd = 
        rdrobust(x = as.numeric(datasets_eff[[i]][, "blmrv"]),
                 y = datasets_eff[[i]][, outcomes[i]],
                 p = polys[j],
                 kernel = kerns[l])
      
      kern_list[l, 1] = outrd$coef[3]
      kern_list[l, 2] = outrd$se[3]
      kern_list[l, 3] = outrd$pv[3]
      kern_list[l, 4] = polys[j]
      kern_list[l, 5] = kerns[l]
      kern_list[l, 6] = outcomes[i]
      
    }
    
    poly_list[[j]] = kern_list %>% 
      as.data.frame %>% 
      `colnames<-` (c("est", "se", "pv", "polys", "kerns", "outcomes")) %>% 
      mutate(est = as.numeric(as.character(est)),
             se = as.numeric(as.character(se)),
             pv = as.numeric(as.character(pv)),
             polys = as.numeric(as.character(polys)),
             kerns = as.character(kerns),
             outcomes = as.character(outcomes))
    
  }
  
  dataset_list[[i]] = poly_list
  
}


dataset_list = dataset_list[1:2]

eff_rdit_res = 
  dataset_list %>% 
  unlist(recursive = FALSE) %>% 
  do.call(rbind.data.frame, .) %>% 
  mutate(outcomes = factor(outcomes,
                           levels = c("arrest",
                                      "hit"),
                           labels = c("A. Arrest Rate | Terry Stop",
                                      "B. Hit Rate | Terry Stop")),
         polys = factor(polys, levels = c(0, 1, 2),
                        labels = c("DIM", "Linear", "Quadratic")),
         kerns = factor(kerns, levels = c("tri", "uni", "epa"),
                        labels = c("Triangular",
                                   "Uniform",
                                   "Epanechnikov")))

eff_rdit_res %>% 
  ggplot() + 
  geom_point(aes(x = polys, y = est, col = kerns),
             position = position_dodge(.4)) +
  geom_errorbar(aes(x = polys, 
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    col = kerns),
                position = position_dodge(.4),
                width = 0,
                size = .2) + 
  geom_errorbar(aes(x = polys, 
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    col = kerns),
                position = position_dodge(.4),
                width = 0,
                size = .4) + 
  geom_hline(yintercept = 0) + 
  scale_color_grey(start = 0, end = .6) + 
  facet_wrap(~ outcomes, scales = "free") + 
  labs(x = "Running Variable Specification",
       y = "Coefficient\n(BLM Protest)",
       col = "Kernel") + 
  theme_tufte()

ggsave(plot = last_plot(), width = 8, height = 2, filename = "resrdit2.png")

#### analyses --- efficiency (closer bandwidth) #### 

# datasets

datasets_eff = list(terry_arr_ts %>% as.data.frame,
                      terry_hit_ts %>% as.data.frame)

# outcomes

outcomes = c("arrest", "hit")

# polynomials

polys = c(0, 1, 2)

# kernels 

kerns = c("tri", "uni", "epa")

# loop for analysis 

dataset_list = as.list(rep(NA, length(datasets_eff)))

for (i in 1:length(datasets_eff)) {
  
  poly_list = as.list(rep(NA, length(polys)))
  
  for (j in 1:length(polys)) {
    
    print(paste0("J Iteration ", j))
    
    kern_list = as.list(rep(NA, length(kerns)))
    
    for (l in 1:length(kerns)) {
      
      print(paste0("L Iteration ", l))
      
      mat_bw = matrix(NA, 100, 7)
      
      for(z in 10:100) {
        
        print(paste0("Z Iteration ", z))
        
        outrd = 
          rdrobust(x = as.numeric(datasets_eff[[i]][, "blmrv"]),## correct rv?
                   y = datasets_eff[[i]][, outcomes[i]],
                   p = polys[j],
                   h = z,
                   kernel = kerns[l])
        
        mat_bw[z, 1] = outrd$coef[3]
        mat_bw[z, 2] = outrd$se[3]
        mat_bw[z, 3] = outrd$pv[3]
        mat_bw[z, 4] = polys[j]
        mat_bw[z, 5] = kerns[l]
        mat_bw[z, 6] = outcomes[i]
        mat_bw[z, 7] = z
        
      }
      
      kern_list[[l]] = mat_bw %>%
        as.data.frame() %>%
        `colnames<-` (c("Coefficient", "SE", "Pvalue", "Polynomial", "Kernel", "Outcome", "BW"))%>% 
        mutate(Coefficient = as.numeric(as.character(Coefficient)),
               SE = as.numeric(as.character(SE)),
               Pvalue = as.numeric(as.character(Pvalue)),
               Polynomial = as.numeric(as.character(Polynomial)),
               Kernel = as.character(Kernel),
               Outcome = as.character(Outcome),
               BW = as.character(BW))
    }
    
    poly_list[[j]] = kern_list
    
  }
  
  dataset_list[[i]] = poly_list
  
}

eff_rdit_res = dataset_list %>% 
  unlist(recursive = FALSE) %>%
  unlist(recursive = FALSE) %>%
  do.call(rbind.data.frame, .) %>% 
  mutate(Outcome = factor(Outcome,
                          levels = c("arrest",
                                     "hit"),
                          labels = c("A. Arrest Rate | Terry Stop",
                                     "B. Hit Rate | Terry Stop")),
         Polynomial = factor(Polynomial, levels = c(0, 1, 2),
                             labels = c("DIM", "Linear", "Quadratic")),
         Kernel = factor(Kernel, levels = c("tri", "uni", "epa"),
                         labels = c("Triangular",
                                    "Uniform",
                                    "Epanechnikov"))) %>% 
  filter(!is.na(BW))

eff_rdit_res = eff_rdit_res %>% 
  mutate(kern_poly = factor(paste0(Kernel, ", ", Polynomial), levels = c("Triangular, DIM", 
                                                                         "Triangular, Linear",
                                                                         "Triangular, Quadratic",
                                                                         "Uniform, DIM",
                                                                         "Uniform, Linear",
                                                                         "Uniform, Quadratic",
                                                                         "Epanechnikov, DIM",
                                                                         "Epanechnikov, Linear",
                                                                         "Epanechnikov, Quadratic")))

eff_rdit_res$BW = as.numeric(eff_rdit_res$BW)

eff_rdit_res %>% filter(Outcome == "B. Hit Rate | Terry Stop") %>% 
  ggplot + geom_point(aes(x = BW, y = Coefficient)) + 
  geom_errorbar(aes(x = BW, ymin = Coefficient - 1.96*SE, ymax = Coefficient + 1.96*SE), width = 0) + 
  facet_wrap(~kern_poly, ncol = 3) + 
  geom_hline(yintercept = 0) + 
  theme_tufte() +
  labs(x = "Bandwidth", y = "RDIT Coefficient (BLM)", title = "Hit rate") +
  scale_x_continuous(breaks = seq(from = 10, to = 100, by = 10))

ggsave(plot = last_plot(), filename = "bw_plot_eff.png", width = 8, height = 6)

qual_rdit_eff = eff_rdit_res %>% filter(Outcome == "B. Hit Rate | Terry Stop") %>%
  filter(Polynomial == "Linear") %>%
  filter(Kernel == "Uniform") %>%
  ggplot + geom_point(aes(x = BW, y = Coefficient)) + 
  geom_errorbar(aes(x = BW, ymin = Coefficient - 1.96*SE, ymax = Coefficient + 1.96*SE), width = 0) + 
  facet_wrap(~kern_poly, ncol = 3) + 
  geom_hline(yintercept = 0) + 
  theme_tufte() +
  labs(x = "Bandwidth", y = "RDIT Coefficient (BLM)", title = "B. Hit rate") +
  scale_x_continuous(breaks = seq(from = 10, to = 100, by = 10))

eff_rdit_res %>% filter(Outcome == "A. Arrest Rate | Terry Stop") %>% 
  ggplot + geom_point(aes(x = BW, y = Coefficient)) + 
  geom_errorbar(aes(x = BW, ymin = Coefficient - 1.96*SE, ymax = Coefficient + 1.96*SE), width = 0) + 
  facet_wrap(~kern_poly, ncol = 3) + 
  geom_hline(yintercept = 0) + 
  theme_tufte() +
  labs(x = "Bandwidth", y = "RDIT Coefficient (BLM)", title = "Arrest rate") +
  scale_x_continuous(breaks = seq(from = 10, to = 100, by = 10))

ggsave(plot = last_plot(), filename = "bw_plot_arrest.png", width = 8, height = 6)

qual_rdit_arrest = eff_rdit_res %>% filter(Outcome == "A. Arrest Rate | Terry Stop") %>%
  filter(Polynomial == "Linear") %>%
  filter(Kernel == "Uniform") %>%
  ggplot + geom_point(aes(x = BW, y = Coefficient)) + 
  geom_errorbar(aes(x = BW, ymin = Coefficient - 1.96*SE, ymax = Coefficient + 1.96*SE), width = 0) + 
  facet_wrap(~kern_poly, ncol = 3) + 
  geom_hline(yintercept = 0) + 
  theme_tufte() +
  labs(x = "Bandwidth", y = "RDIT Coefficient (BLM)", title = "A. Arrest rate") +
  scale_x_continuous(breaks = seq(from = 10, to = 100, by = 10))

#### analyses --- disparities #### 

# datasets

terry_arr_ts = 
  terry %>% 
  mutate(count = 1,
         r_unk = ifelse(subject_race == "Unknown" | subject_race == "-", 1, 0),
         r_wht = ifelse(subject_race == "White", 1, 0),
         r_asn = ifelse(subject_race == "Asian", 1, 0),
         r_blk = ifelse(subject_race == "Black or African American", 1, 0)) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            arrest = mean(arrest, na.rm = TRUE),
            r_wht = sum(r_wht, na.rm = TRUE),
            r_blk = sum(r_blk, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2019-06-01")) %>% 
  mutate(wht_rate = (r_wht / (.673 * 737015)) * 10000,
         blk_rate = (r_blk / (.073 * 737015)) * 10000) %>% 
  mutate(rr_bw = blk_rate / wht_rate) %>% 
  mutate(rr_bw = ifelse(rr_bw == Inf | is.na(rr_bw), NA, rr_bw))

# running variable

terry_arr_ts$blmrv = 
  terry_arr_ts$date - min(terry_arr_ts$date[terry_arr_ts$blm == 1])

# datasets

datasets_eff = list(terry_arr_ts %>% as.data.frame)

# outcomes

outcomes = c("rr_bw")

# polynomials

polys = c(0, 1, 2)

# kernels 

kerns = c("tri", "uni", "epa")

# loop for analysis 

dataset_list = as.list(rep(NA, length(datasets_eff)))

for (i in 1:length(datasets_eff)) {
  
  poly_list = as.list(rep(NA, length(polys)))
  
  for (j in 1:length(polys)) {
    
    print(paste0("J Iteration ", j))
    
    kern_list = matrix(nrow = length(kerns), ncol = 6)
    
    for (l in 1:length(kerns)) {
      
      print(paste0("L Iteration ", l))
      
      outrd = 
        rdrobust(x = as.numeric(datasets_eff[[i]][, "blmrv"]),
                 y = datasets_eff[[i]][, outcomes[i]],
                 p = polys[j],
                 kernel = kerns[l])
      
      kern_list[l, 1] = outrd$coef[3]
      kern_list[l, 2] = outrd$se[3]
      kern_list[l, 3] = outrd$pv[3]
      kern_list[l, 4] = polys[j]
      kern_list[l, 5] = kerns[l]
      kern_list[l, 6] = outcomes[i]
      
    }
    
    poly_list[[j]] = kern_list %>% 
      as.data.frame %>% 
      `colnames<-` (c("est", "se", "pv", "polys", "kerns", "outcomes")) %>% 
      mutate(est = as.numeric(as.character(est)),
             se = as.numeric(as.character(se)),
             pv = as.numeric(as.character(pv)),
             polys = as.numeric(as.character(polys)),
             kerns = as.character(kerns),
             outcomes = as.character(outcomes))
    
  }
  
  dataset_list[[i]] = poly_list
  
}

eff_rdit_res = 
  dataset_list[[1]] %>% 
  do.call(rbind.data.frame, .) %>% 
  mutate(outcomes = factor(outcomes,
                           levels = c("rr_bw"),
                           labels = c("Black/White Rate Ratio")),
         polys = factor(polys, levels = c(0, 1, 2),
                        labels = c("DIM", "Linear", "Quadratic")),
         kerns = factor(kerns, levels = c("tri", "uni", "epa"),
                        labels = c("Triangular",
                                   "Uniform",
                                   "Epanechnikov")))

eff_rdit_res %>% 
  ggplot() + 
  geom_point(aes(x = polys, y = est, col = kerns),
             position = position_dodge(.4)) +
  geom_errorbar(aes(x = polys, 
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    col = kerns),
                position = position_dodge(.4),
                width = 0,
                size = .2) + 
  geom_errorbar(aes(x = polys, 
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    col = kerns),
                position = position_dodge(.4),
                width = 0,
                size = .4) + 
  geom_hline(yintercept = 0) + 
  scale_color_grey(start = 0, end = .6) + 
  facet_wrap(~ outcomes, scales = "free") + 
  labs(x = "Running Variable Specification",
       y = "Coefficient\n(BLM Protest)",
       col = "Kernel") + 
  theme_tufte()

ggsave(plot = last_plot(), width = 4, height = 2, filename = "resdisp.png")

#### analyses --- disparities (closer bandwidth) ####

# datasets

terry_arr_ts = 
  terry %>% 
  mutate(count = 1,
         r_unk = ifelse(subject_race == "Unknown" | subject_race == "-", 1, 0),
         r_wht = ifelse(subject_race == "White", 1, 0),
         r_asn = ifelse(subject_race == "Asian", 1, 0),
         r_blk = ifelse(subject_race == "Black or African American", 1, 0)) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            arrest = mean(arrest, na.rm = TRUE),
            r_wht = sum(r_wht, na.rm = TRUE),
            r_blk = sum(r_blk, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2019-06-01")) %>% 
  mutate(wht_rate = (r_wht / (.673 * 737015)) * 10000,
         blk_rate = (r_blk / (.073 * 737015)) * 10000) %>% 
  mutate(rr_bw = blk_rate / wht_rate) %>% 
  mutate(rr_bw = ifelse(rr_bw == Inf | is.na(rr_bw), NA, rr_bw))

# running variable

terry_arr_ts$blmrv = 
  terry_arr_ts$date - min(terry_arr_ts$date[terry_arr_ts$blm == 1])

# datasets

datasets_eff = list(terry_arr_ts %>% as.data.frame)

# outcomes

outcomes = c("rr_bw")

# polynomials

polys = c(0, 1, 2)

# kernels 

kerns = c("tri", "uni", "epa")

# loop for analysis 

dataset_list = as.list(rep(NA, length(datasets_eff)))

for (i in 1:length(datasets_eff)) {
  
  poly_list = as.list(rep(NA, length(polys)))
  
  for (j in 1:length(polys)) {
    
    print(paste0("J Iteration ", j))
    
    kern_list = as.list(rep(NA, length(kerns)))
    
    for (l in 1:length(kerns)) {
      
      print(paste0("L Iteration ", l))
      
      mat_bw = matrix(NA, 100, 7)
      
      for(z in 10:100) {
        
        print(paste0("Z Iteration ", z))
        
        outrd = 
          rdrobust(x = as.numeric(datasets_eff[[i]][, "blmrv"]),## correct rv?
                   y = datasets_eff[[i]][, outcomes[i]],
                   p = polys[j],
                   h = z,
                   kernel = kerns[l])
        
        mat_bw[z, 1] = outrd$coef[3]
        mat_bw[z, 2] = outrd$se[3]
        mat_bw[z, 3] = outrd$pv[3]
        mat_bw[z, 4] = polys[j]
        mat_bw[z, 5] = kerns[l]
        mat_bw[z, 6] = outcomes[i]
        mat_bw[z, 7] = z
        
      }
      
      kern_list[[l]] = mat_bw %>%
        as.data.frame() %>%
        `colnames<-` (c("Coefficient", "SE", "Pvalue", "Polynomial", "Kernel", "Outcome", "BW"))%>% 
        mutate(Coefficient = as.numeric(as.character(Coefficient)),
               SE = as.numeric(as.character(SE)),
               Pvalue = as.numeric(as.character(Pvalue)),
               Polynomial = as.numeric(as.character(Polynomial)),
               Kernel = as.character(Kernel),
               Outcome = as.character(Outcome),
               BW = as.character(BW))
    }
    
    poly_list[[j]] = kern_list
    
  }
  
  dataset_list[[i]] = poly_list
  
}

eff_rdit_res = dataset_list %>% 
  unlist(recursive = FALSE) %>%
  unlist(recursive = FALSE) %>%
  do.call(rbind.data.frame, .) %>% 
  mutate(Outcome = factor(Outcome,
                          levels = c("rr_bw"),
                          labels = c("Black/White Rate Ratio")),
         Polynomial = factor(Polynomial, levels = c(0, 1, 2),
                             labels = c("DIM", "Linear", "Quadratic")),
         Kernel = factor(Kernel, levels = c("tri", "uni", "epa"),
                         labels = c("Triangular",
                                    "Uniform",
                                    "Epanechnikov"))) %>% 
  filter(!is.na(BW))

eff_rdit_res = eff_rdit_res %>% 
  mutate(kern_poly = factor(paste0(Kernel, ", ", Polynomial), levels = c("Triangular, DIM", 
                                                                         "Triangular, Linear",
                                                                         "Triangular, Quadratic",
                                                                         "Uniform, DIM",
                                                                         "Uniform, Linear",
                                                                         "Uniform, Quadratic",
                                                                         "Epanechnikov, DIM",
                                                                         "Epanechnikov, Linear",
                                                                         "Epanechnikov, Quadratic")))

eff_rdit_res$BW = as.numeric(eff_rdit_res$BW)

eff_rdit_res %>% filter(Outcome == "Black/White Rate Ratio") %>% 
  ggplot + geom_point(aes(x = BW, y = Coefficient)) + 
  geom_errorbar(aes(x = BW, ymin = Coefficient - 1.96*SE, ymax = Coefficient + 1.96*SE), width = 0) + 
  facet_wrap(~kern_poly, ncol = 3) + 
  geom_hline(yintercept = 0) + 
  theme_tufte() +
  labs(x = "Bandwidth", y = "RDIT Coefficient (BLM)", title = "Black/White disparity") +
  scale_x_continuous(breaks = seq(from = 10, to = 100, by = 10))

ggsave(plot = last_plot(), filename = "bw_plot_disparity.png", width = 8, height = 6)

qual_rdit_ratio = eff_rdit_res %>% filter(Outcome == "Black/White Rate Ratio") %>%
  filter(Polynomial == "Linear") %>%
  filter(Kernel == "Uniform") %>%
  ggplot + geom_point(aes(x = BW, y = Coefficient)) + 
  geom_errorbar(aes(x = BW, ymin = Coefficient - 1.96*SE, ymax = Coefficient + 1.96*SE), width = 0) + 
  facet_wrap(~kern_poly, ncol = 3) + 
  geom_hline(yintercept = 0) + 
  theme_tufte() +
  labs(x = "Bandwidth", y = "RDIT Coefficient (BLM)", title = "C. Black/White disparity") +
  scale_x_continuous(breaks = seq(from = 10, to = 100, by = 10))

qual_uni_lin = arrangeGrob(qual_rdit_arrest, qual_rdit_eff, qual_rdit_ratio, depol_plot3, ncol = 2)
ggsave(filename = "qual_uni_lin.png", plot = qual_uni_lin, width = 8, height = 4)

#### analyses --- crime ####

datasets_crm = list(crime_ts %>% as.data.frame %>% filter(date >= as.Date("2016-01-01")),
                    crime_ts %>% as.data.frame %>% filter(date >= as.Date("2016-01-01")),
                    crime_ts %>% as.data.frame %>% filter(date >= as.Date("2016-01-01")),
                    crime_ts %>% as.data.frame %>% filter(date >= as.Date("2016-01-01")))

# tab for against person crime 

prop.table(table(crime$offense_parent_group[crime$date >= as.Date("2016-01-01") & crime$cat_person == 1]))[
  order(prop.table(table(crime$offense_parent_group[crime$date >= as.Date("2016-01-01") & crime$cat_person == 1])))
]

# tab for against property crime 

prop.table(table(crime$offense_parent_group[crime$date >= as.Date("2016-01-01") & crime$cat_property == 1]))[
  order(prop.table(table(crime$offense_parent_group[crime$date >= as.Date("2016-01-01") & crime$cat_property == 1])))
  ]

# tab for against society crime 

prop.table(table(crime$offense_parent_group[crime$date >= as.Date("2016-01-01") & crime$cat_society == 1]))[
  order(prop.table(table(crime$offense_parent_group[crime$date >= as.Date("2016-01-01") & crime$cat_society == 1])))
  ]

# outcomes

outcomes = c("count", "cat_property", "cat_society", "cat_person")

# polynomials

polys = c(0, 1, 2)

# kernels 

kerns = c("tri", "uni", "epa")

# loop for analysis 

dataset_list = as.list(rep(NA, length(datasets_crm)))

for (i in 1:length(datasets_crm)) {
  
  poly_list = as.list(rep(NA, length(polys)))
  
  for (j in 1:length(polys)) {
    
    print(paste0("J Iteration ", j))
    
    kern_list = matrix(nrow = length(kerns), ncol = 6)
    
    for (l in 1:length(kerns)) {
      
      print(paste0("L Iteration ", l))
      
      outrd = 
        rdrobust(x = as.numeric(datasets_crm[[i]][, "blmrv"]),
                 y = datasets_crm[[i]][, outcomes[i]],
                 p = polys[j],
                 kernel = kerns[l])
      
      kern_list[l, 1] = outrd$coef[3]
      kern_list[l, 2] = outrd$se[3]
      kern_list[l, 3] = outrd$pv[3]
      kern_list[l, 4] = polys[j]
      kern_list[l, 5] = kerns[l]
      kern_list[l, 6] = outcomes[i]
      
    }
    
    poly_list[[j]] = kern_list %>% 
      as.data.frame %>% 
      `colnames<-` (c("est", "se", "pv", "polys", "kerns", "outcomes")) %>% 
      mutate(est = as.numeric(as.character(est)),
             se = as.numeric(as.character(se)),
             pv = as.numeric(as.character(pv)),
             polys = as.numeric(as.character(polys)),
             kerns = as.character(kerns),
             outcomes = as.character(outcomes))
    
  }
  
  dataset_list[[i]] = poly_list
  
}

crm_rdit_res = 
  dataset_list %>% 
  unlist(recursive = FALSE) %>% 
  do.call(rbind.data.frame, .) %>% 
  mutate(outcomes = factor(outcomes,
                           levels = c("count",
                                      "cat_property",
                                      "cat_society",
                                      "cat_person"),
                           labels = c("A. All Crimes",
                                      "B. Property Crimes",
                                      "C. Against Society Crimes",
                                      "D. Against Person Crimes")),
         polys = factor(polys, levels = c(0, 1, 2),
                        labels = c("DIM", "Linear", "Quadratic")),
         kerns = factor(kerns, levels = c("tri", "uni", "epa"),
                        labels = c("Triangular",
                                   "Uniform",
                                   "Epanechnikov")))

crm_rdit_res %>% 
  ggplot() + 
  geom_point(aes(x = polys, y = est, col = kerns),
             position = position_dodge(.6)) +
  geom_errorbar(aes(x = polys, 
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    col = kerns),
                position = position_dodge(.6),
                width = 0,
                size = .2) + 
  geom_errorbar(aes(x = polys, 
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    col = kerns),
                position = position_dodge(.6),
                width = 0,
                size = .4) + 
  geom_hline(yintercept = 0) + 
  scale_color_grey(start = 0, end = .6) + 
  facet_wrap(~ outcomes, scales = "free", ncol = 4) + 
  labs(x = "Running Variable Specification",
       y = "Coefficient\n(BLM Protest)",
       col = "Kernel") + 
  theme_tufte(base_size = 9)

ggsave(plot = last_plot(), width = 8, height = 2, filename = "resrdit3.png")

#### analyses ---- crime (closer bandwidth) ####

datasets_crm = list(crime_ts %>% as.data.frame %>% filter(date >= as.Date("2016-01-01")),
                    crime_ts %>% as.data.frame %>% filter(date >= as.Date("2016-01-01")),
                    crime_ts %>% as.data.frame %>% filter(date >= as.Date("2016-01-01")),
                    crime_ts %>% as.data.frame %>% filter(date >= as.Date("2016-01-01")))

# tab for against person crime 

prop.table(table(crime$offense_parent_group[crime$date >= as.Date("2016-01-01") & crime$cat_person == 1]))[
  order(prop.table(table(crime$offense_parent_group[crime$date >= as.Date("2016-01-01") & crime$cat_person == 1])))
]

# tab for against property crime 

prop.table(table(crime$offense_parent_group[crime$date >= as.Date("2016-01-01") & crime$cat_property == 1]))[
  order(prop.table(table(crime$offense_parent_group[crime$date >= as.Date("2016-01-01") & crime$cat_property == 1])))
]

# tab for against society crime 

prop.table(table(crime$offense_parent_group[crime$date >= as.Date("2016-01-01") & crime$cat_society == 1]))[
  order(prop.table(table(crime$offense_parent_group[crime$date >= as.Date("2016-01-01") & crime$cat_society == 1])))
]

# outcomes

outcomes = c("count", "cat_property", "cat_society", "cat_person")

# polynomials

polys = c(0, 1, 2)

# kernels 

kerns = c("tri", "uni", "epa")

# loop for analysis 

dataset_list = as.list(rep(NA, length(datasets_crm)))

for (i in 1:length(datasets_crm)) {
  
  poly_list = as.list(rep(NA, length(polys)))
  
  for (j in 1:length(polys)) {
    
    print(paste0("J Iteration ", j))
    
    kern_list = as.list(rep(NA, length(kerns)))
    
    for (l in 1:length(kerns)) {
      
      print(paste0("L Iteration ", l))
      
      mat_bw = matrix(NA, 100, 7)
      
      for(z in 10:100) {
        
        print(paste0("Z Iteration ", z))
        
        outrd = 
          rdrobust(x = as.numeric(datasets_crm[[i]][, "blmrv"]),## correct rv?
                   y = datasets_crm[[i]][, outcomes[i]],
                   p = polys[j],
                   h = z,
                   kernel = kerns[l])
        
        mat_bw[z, 1] = outrd$coef[3]
        mat_bw[z, 2] = outrd$se[3]
        mat_bw[z, 3] = outrd$pv[3]
        mat_bw[z, 4] = polys[j]
        mat_bw[z, 5] = kerns[l]
        mat_bw[z, 6] = outcomes[i]
        mat_bw[z, 7] = z
        
      }
      
      kern_list[[l]] = mat_bw %>%
        as.data.frame() %>%
        `colnames<-` (c("Coefficient", "SE", "Pvalue", "Polynomial", "Kernel", "Outcome", "BW"))%>% 
        mutate(Coefficient = as.numeric(as.character(Coefficient)),
               SE = as.numeric(as.character(SE)),
               Pvalue = as.numeric(as.character(Pvalue)),
               Polynomial = as.numeric(as.character(Polynomial)),
               Kernel = as.character(Kernel),
               Outcome = as.character(Outcome),
               BW = as.character(BW))
    }
    
    poly_list[[j]] = kern_list
    
  }
  
  dataset_list[[i]] = poly_list
  
}

crm_rdit_res = dataset_list %>% 
  unlist(recursive = FALSE) %>%
  unlist(recursive = FALSE) %>%
  do.call(rbind.data.frame, .) %>% 
  mutate(Outcome = factor(Outcome,
                          levels = c("count",
                                     "cat_property",
                                     "cat_society",
                                     "cat_person"),
                          labels = c("A. All Crimes",
                                     "B. Property Crimes",
                                     "C. Against Society Crimes",
                                     "D. Against Person Crimes")),
         Polynomial = factor(Polynomial, levels = c(0, 1, 2),
                             labels = c("DIM", "Linear", "Quadratic")),
         Kernel = factor(Kernel, levels = c("tri", "uni", "epa"),
                         labels = c("Triangular",
                                    "Uniform",
                                    "Epanechnikov"))) %>% 
  filter(!is.na(BW))

crm_rdit_res = crm_rdit_res %>% 
  mutate(kern_poly = factor(paste0(Kernel, ", ", Polynomial), levels = c("Triangular, DIM", 
                                                                         "Triangular, Linear",
                                                                         "Triangular, Quadratic",
                                                                         "Uniform, DIM",
                                                                         "Uniform, Linear",
                                                                         "Uniform, Quadratic",
                                                                         "Epanechnikov, DIM",
                                                                         "Epanechnikov, Linear",
                                                                         "Epanechnikov, Quadratic")))

crm_rdit_res$BW = as.numeric(crm_rdit_res$BW)

crm_rdit_res %>% filter(Outcome == "A. All Crimes") %>% 
  ggplot + geom_point(aes(x = BW, y = Coefficient)) + 
  geom_errorbar(aes(x = BW, ymin = Coefficient - 1.96*SE, ymax = Coefficient + 1.96*SE), width = 0) + 
  facet_wrap(~kern_poly, ncol = 3) + 
  geom_hline(yintercept = 0) + 
  theme_tufte() +
  labs(x = "Bandwidth", y = "RDIT Coefficient (BLM)", title = "Crime rate (all crimes)") +
  scale_x_continuous(breaks = seq(from = 10, to = 100, by = 10))

ggsave(plot = last_plot(), filename = "bw_plot_allcrime.png", width = 8, height = 6)

crm_rdit_res %>% filter(Outcome == "B. Property Crimes") %>% 
  ggplot + geom_point(aes(x = BW, y = Coefficient)) + 
  geom_errorbar(aes(x = BW, ymin = Coefficient - 1.96*SE, ymax = Coefficient + 1.96*SE), width = 0) + 
  facet_wrap(~kern_poly, ncol = 3) + 
  geom_hline(yintercept = 0) + 
  theme_tufte() +
  labs(x = "Bandwidth", y = "RDIT Coefficient (BLM)", title = "Crime rate (property crimes)") +
  scale_x_continuous(breaks = seq(from = 10, to = 100, by = 10))

ggsave(plot = last_plot(), filename = "bw_plot_propcrime.png", width = 8, height = 6)

crm_rdit_res %>% filter(Outcome == "C. Against Society Crimes") %>% 
  ggplot + geom_point(aes(x = BW, y = Coefficient)) + 
  geom_errorbar(aes(x = BW, ymin = Coefficient - 1.96*SE, ymax = Coefficient + 1.96*SE), width = 0) + 
  facet_wrap(~kern_poly, ncol = 3) + 
  geom_hline(yintercept = 0) + 
  theme_tufte() +
  labs(x = "Bandwidth", y = "RDIT Coefficient (BLM)", title = "Crime rate (against society crimes)") +
  scale_x_continuous(breaks = seq(from = 10, to = 100, by = 10))

ggsave(plot = last_plot(), filename = "bw_plot_soccrime.png", width = 8, height = 6)

crm_rdit_res %>% filter(Outcome == "D. Against Person Crimes") %>% 
  ggplot + geom_point(aes(x = BW, y = Coefficient)) + 
  geom_errorbar(aes(x = BW, ymin = Coefficient - 1.96*SE, ymax = Coefficient + 1.96*SE), width = 0) + 
  facet_wrap(~kern_poly, ncol = 3) + 
  geom_hline(yintercept = 0) + 
  theme_tufte() +
  labs(x = "Bandwidth", y = "RDIT Coefficient (BLM)", title = "Crime rate (against person crimes)") +
  scale_x_continuous(breaks = seq(from = 10, to = 100, by = 10))

ggsave(plot = last_plot(), filename = "bw_plot_personcrime.png", width = 8, height = 6)

#### Crime -- plots for Yale paper ####

crm_plot_all = crm_rdit_res %>% filter(Outcome == "A. All Crimes") %>%
  filter(Polynomial == "Linear") %>%
  filter(Kernel == "Uniform") %>%
  ggplot + geom_point(aes(x = BW, y = Coefficient)) + 
  geom_errorbar(aes(x = BW, ymin = Coefficient - 1.96*SE, ymax = Coefficient + 1.96*SE), width = 0) + 
  facet_wrap(~kern_poly, ncol = 3) + 
  geom_hline(yintercept = 0) + 
  theme_tufte() +
  labs(x = "Bandwidth", y = "RDIT Coefficient (BLM)", title = "A. All Crimes") +
  scale_x_continuous(breaks = seq(from = 10, to = 100, by = 10))

crm_plot_prop = crm_rdit_res %>% filter(Outcome == "B. Property Crimes") %>%
  filter(Polynomial == "Linear") %>%
  filter(Kernel == "Uniform") %>%
  ggplot + geom_point(aes(x = BW, y = Coefficient)) + 
  geom_errorbar(aes(x = BW, ymin = Coefficient - 1.96*SE, ymax = Coefficient + 1.96*SE), width = 0) + 
  facet_wrap(~kern_poly, ncol = 3) + 
  geom_hline(yintercept = 0) + 
  theme_tufte() +
  labs(x = "Bandwidth", y = "RDIT Coefficient (BLM)", title = "C. Against Property Crimes") +
  scale_x_continuous(breaks = seq(from = 10, to = 100, by = 10))

crm_plot_soc = crm_rdit_res %>% filter(Outcome == "C. Against Society Crimes") %>%
  filter(Polynomial == "Linear") %>%
  filter(Kernel == "Uniform") %>%
  ggplot + geom_point(aes(x = BW, y = Coefficient)) + 
  geom_errorbar(aes(x = BW, ymin = Coefficient - 1.96*SE, ymax = Coefficient + 1.96*SE), width = 0) + 
  facet_wrap(~kern_poly, ncol = 3) + 
  geom_hline(yintercept = 0) + 
  theme_tufte() +
  labs(x = "Bandwidth", y = "RDIT Coefficient (BLM)", title = "D. Against Society Crimes") +
  scale_x_continuous(breaks = seq(from = 10, to = 100, by = 10))

crm_plot_pers = crm_rdit_res %>% filter(Outcome == "D. Against Person Crimes") %>%
  filter(Polynomial == "Linear") %>%
  filter(Kernel == "Uniform") %>%
  ggplot + geom_point(aes(x = BW, y = Coefficient)) + 
  geom_errorbar(aes(x = BW, ymin = Coefficient - 1.96*SE, ymax = Coefficient + 1.96*SE), width = 0) + 
  facet_wrap(~kern_poly, ncol = 3) + 
  geom_hline(yintercept = 0) + 
  theme_tufte() +
  labs(x = "Bandwidth", y = "RDIT Coefficient (BLM)", title = "B. Against Person Crimes") +
  scale_x_continuous(breaks = seq(from = 10, to = 100, by = 10))

crm_grob_paper = arrangeGrob(crm_plot_all, crm_plot_pers, crm_plot_prop, crm_plot_soc, ncol = 2)
ggsave(filename = "crm_uni_lin.png", plot = crm_grob_paper, width = 8, height = 4)

#### analyses, close to disc --- depol #### 

#### analyses --- depolicing --- semi-donut hole (to demonstrate effect persistence) #### 

# datasets

datasets_depol = list(terry_ts %>% as.data.frame %>% mutate(count_terry = count) %>% filter(date >= as.Date("2016-01-01")),
                      call_ts %>% as.data.frame %>% mutate(count_calls = count) %>% filter(date >= as.Date("2016-01-01")),
                      call_ts %>% as.data.frame %>% filter(date >= as.Date("2016-01-01")),
                      call_ts2 %>% as.data.frame %>% mutate(count_calls_civ = count) %>% filter(date >= as.Date("2016-01-01")),
                      call_ts3 %>% as.data.frame %>% mutate(count_calls_off = count) %>% filter(date >= as.Date("2016-01-01")))

# outcomes

outcomes = c("count_terry", "count_calls", "response_time", "count_calls_civ", "count_calls_off")

# polynomials

polys = c(0, 1, 2)

# kernels 

kerns = c("tri", "uni", "epa")

# loop for analysis 

dataset_list = as.list(rep(NA, length(datasets_depol)))

to_cut_rhs = 1:317

for (i in 1:length(datasets_depol)) {
  
  print(paste0("I Iteration ", i))
  
  poly_list = as.list(rep(NA, length(polys)))
  
  for (j in 1:length(polys)) {
    
    print(paste0("J Iteration ", j))
    
    kern_list = as.list(rep(NA, length(kerns)))
    
    for (l in 1:length(kerns)) {
      
      print(paste0("L Iteration ", l))
      
      dnut_mat = matrix(nrow = length(to_cut_rhs), ncol = 7)
      
      for (z in 1:length(to_cut_rhs)) {
        
        print(paste0("Z Iteration ", z))
        
        
        cut_data = 
          datasets_depol[[i]] %>% 
          filter(blmrv < 0 | blmrv >= z) %>% 
          mutate(blmrv = ifelse(blmrv > 0, blmrv - z, blmrv))
        
        outrd = 
          rdrobust(x = as.numeric(cut_data[, "blmrv"]),
                   y = cut_data[, outcomes[i]],
                   p = polys[j],
                   kernel = kerns[l],
                   h = 50)
        
        dnut_mat[z, 1] = outrd$coef[3]
        dnut_mat[z, 2] = outrd$se[3]
        dnut_mat[z, 3] = outrd$pv[3]
        dnut_mat[z, 4] = polys[j]
        dnut_mat[z, 5] = kerns[l]
        dnut_mat[z, 6] = outcomes[i]
        dnut_mat[z, 7] = to_cut_rhs[z]
        
      }
      
      kern_list[[l]] = dnut_mat %>% 
        as.data.frame %>% 
        `colnames<-` (c("est", "se", "pv", "polys", "kerns", "outcomes", "cut")) %>% 
        mutate(est = as.numeric(as.character(est)),
               se = as.numeric(as.character(se)),
               pv = as.numeric(as.character(pv)),
               polys = as.numeric(as.character(polys)),
               kerns = as.character(kerns),
               outcomes = as.character(outcomes),
               cut = as.numeric(as.character(cut)))
      
    }
    
    poly_list[[j]] = kern_list 
    
  }
  
  dataset_list[[i]] = poly_list
  
}

plot_cut_bw1 = dataset_list %>% 
  unlist(recursive = FALSE) %>% 
  unlist(recursive = FALSE) %>% 
  do.call(rbind.data.frame, .) %>% 
  filter(outcomes == "count_terry") %>% 
  mutate(polys = factor(polys, levels = c(0, 1, 2), labels = c("DIM",
                                                       "Linear",
                                                       "Quadratic")),
         kerns = factor(kerns, levels = c('tri', "uni", "epa"),
                        labels = c("Triangular", "Uniform", "Epanechnikov"))) %>% 
  mutate(kern_poly = paste0(kerns, ", " , polys)) %>% 
  mutate(kern_poly = factor(kern_poly, 
                            levels = c("Triangular, DIM",
                                       "Triangular, Linear",
                                       "Triangular, Quadratic",
                                       "Uniform, DIM",
                                       "Uniform, Linear",
                                       "Uniform, Quadratic",
                                       "Epanechnikov, DIM",
                                       "Epanechnikov, Linear",
                                       "Epanechnikov, Quadratic")))

plot_cut_bw2 = dataset_list %>% 
  unlist(recursive = FALSE) %>% 
  unlist(recursive = FALSE) %>% 
  do.call(rbind.data.frame, .) %>% 
  filter(outcomes == "count_calls") %>% 
  mutate(polys = factor(polys, levels = c(0, 1, 2), labels = c("DIM",
                                                               "Linear",
                                                               "Quadratic")),
         kerns = factor(kerns, levels = c('tri', "uni", "epa"),
                        labels = c("Triangular", "Uniform", "Epanechnikov"))) %>% 
  mutate(kern_poly = paste0(kerns, ", " , polys)) %>% 
  mutate(kern_poly = factor(kern_poly, 
                            levels = c("Triangular, DIM",
                                       "Triangular, Linear",
                                       "Triangular, Quadratic",
                                       "Uniform, DIM",
                                       "Uniform, Linear",
                                       "Uniform, Quadratic",
                                       "Epanechnikov, DIM",
                                       "Epanechnikov, Linear",
                                       "Epanechnikov, Quadratic")))

plot_cut_bw3 = dataset_list %>% 
  unlist(recursive = FALSE) %>% 
  unlist(recursive = FALSE) %>% 
  do.call(rbind.data.frame, .) %>% 
  filter(outcomes == "response_time") %>% 
  mutate(polys = factor(polys, levels = c(0, 1, 2), labels = c("DIM",
                                                               "Linear",
                                                               "Quadratic")),
         kerns = factor(kerns, levels = c('tri', "uni", "epa"),
                        labels = c("Triangular", "Uniform", "Epanechnikov"))) %>% 
  mutate(kern_poly = paste0(kerns, ", " , polys)) %>% 
  mutate(kern_poly = factor(kern_poly, 
                            levels = c("Triangular, DIM",
                                       "Triangular, Linear",
                                       "Triangular, Quadratic",
                                       "Uniform, DIM",
                                       "Uniform, Linear",
                                       "Uniform, Quadratic",
                                       "Epanechnikov, DIM",
                                       "Epanechnikov, Linear",
                                       "Epanechnikov, Quadratic")))

plot_cut_bw4 = dataset_list %>% 
  unlist(recursive = FALSE) %>% 
  unlist(recursive = FALSE) %>% 
  do.call(rbind.data.frame, .) %>% 
  filter(outcomes == "count_calls_off") %>% 
  mutate(polys = factor(polys, levels = c(0, 1, 2), labels = c("DIM",
                                                               "Linear",
                                                               "Quadratic")),
         kerns = factor(kerns, levels = c('tri', "uni", "epa"),
                        labels = c("Triangular", "Uniform", "Epanechnikov"))) %>% 
  mutate(kern_poly = paste0(kerns, ", " , polys)) %>% 
  mutate(kern_poly = factor(kern_poly, 
                            levels = c("Triangular, DIM",
                                       "Triangular, Linear",
                                       "Triangular, Quadratic",
                                       "Uniform, DIM",
                                       "Uniform, Linear",
                                       "Uniform, Quadratic",
                                       "Epanechnikov, DIM",
                                       "Epanechnikov, Linear",
                                       "Epanechnikov, Quadratic")))

plot_cut_bw5 = dataset_list %>% 
  unlist(recursive = FALSE) %>% 
  unlist(recursive = FALSE) %>% 
  do.call(rbind.data.frame, .) %>% 
  filter(outcomes == "count_calls_civ") %>% 
  mutate(polys = factor(polys, levels = c(0, 1, 2), labels = c("DIM",
                                                               "Linear",
                                                               "Quadratic")),
         kerns = factor(kerns, levels = c('tri', "uni", "epa"),
                        labels = c("Triangular", "Uniform", "Epanechnikov"))) %>% 
  mutate(kern_poly = paste0(kerns, ", " , polys)) %>% 
  mutate(kern_poly = factor(kern_poly, 
                            levels = c("Triangular, DIM",
                                       "Triangular, Linear",
                                       "Triangular, Quadratic",
                                       "Uniform, DIM",
                                       "Uniform, Linear",
                                       "Uniform, Quadratic",
                                       "Epanechnikov, DIM",
                                       "Epanechnikov, Linear",
                                       "Epanechnikov, Quadratic")))


plot_cut_bw1 %>% 
  ggplot() + 
  geom_line(aes(x = cut, y = est),
            linetype = 1,
            size = .2) +
  geom_line(aes(x = cut, y = est - 1.96 * se),
            linetype = 2,
            size = .2) +
  geom_line(aes(x = cut, y = est + 1.96 * se),
            linetype = 2,
            size = .2) + 
  geom_hline(yintercept = 0) + 
  labs(x = "Cut Days",
       y = "RDiT Coefficient (BLM)\n(50 Day Bandwidth)",
       title = "Terry Stops") + 
  facet_wrap(~ kern_poly, ncol = 3) + 
  theme_tufte()

ggsave(plot = last_plot(), width = 8, height = 6, filename = "cutdepol1.png")

plot_cut_bw2 %>% 
  ggplot() + 
  geom_line(aes(x = cut, y = est),
            linetype = 1,
            size = .2) +
  geom_line(aes(x = cut, y = est - 1.96 * se),
            linetype = 2,
            size = .2) +
  geom_line(aes(x = cut, y = est + 1.96 * se),
            linetype = 2,
            size = .2) + 
  geom_hline(yintercept = 0) + 
  labs(x = "Cut Days",
       y = "RDiT Coefficient (BLM)\n(50 Day Bandwidth)",
       title = "Total 911 Calls") + 
  facet_wrap(~ kern_poly, ncol = 3) + 
  theme_tufte()

ggsave(plot = last_plot(), width = 8, height = 6, filename = "cutdepol2.png")

plot_cut_bw3 %>% 
  ggplot() + 
  geom_line(aes(x = cut, y = est),
            linetype = 1,
            size = .2) +
  geom_line(aes(x = cut, y = est - 1.96 * se),
            linetype = 2,
            size = .2) +
  geom_line(aes(x = cut, y = est + 1.96 * se),
            linetype = 2,
            size = .2) + 
  geom_hline(yintercept = 0) + 
  labs(x = "Cut Days",
       y = "RDiT Coefficient (BLM)\n(50 Day Bandwidth)",
       title = "Response Times") + 
  facet_wrap(~ kern_poly, ncol = 3) + 
  theme_tufte()

ggsave(plot = last_plot(), width = 8, height = 6, filename = "cutdepol3.png")

plot_cut_bw4 %>% 
  ggplot() + 
  geom_line(aes(x = cut, y = est),
            linetype = 1,
            size = .2) +
  geom_line(aes(x = cut, y = est - 1.96 * se),
            linetype = 2,
            size = .2) +
  geom_line(aes(x = cut, y = est + 1.96 * se),
            linetype = 2,
            size = .2) + 
  geom_hline(yintercept = 0) + 
  labs(x = "Cut Days",
       y = "RDiT Coefficient (BLM)\n(50 Day Bandwidth)",
       title = "Officer Initiated 911 Calls") + 
  facet_wrap(~ kern_poly, ncol = 3) + 
  theme_tufte()

ggsave(plot = last_plot(), width = 8, height = 6, filename = "cutdepol4.png")

plot_cut_bw5 %>% 
  ggplot() + 
  geom_line(aes(x = cut, y = est),
            linetype = 1,
            size = .2) +
  geom_line(aes(x = cut, y = est - 1.96 * se),
            linetype = 2,
            size = .2) +
  geom_line(aes(x = cut, y = est + 1.96 * se),
            linetype = 2,
            size = .2) + 
  geom_hline(yintercept = 0) + 
  labs(x = "Cut Days",
       y = "RDiT Coefficient (BLM)\n(50 Day Bandwidth)",
       title = "Civilian Initiated 911 Calls") + 
  facet_wrap(~ kern_poly, ncol = 3) + 
  theme_tufte()

ggsave(plot = last_plot(), width = 8, height = 6, filename = "cutdepol5.png")

#### geo-analysis ####

terry_geo = 
  terry_merged_use %>%
  group_by(beat) %>%
  summarise(p_nonwhite = mean(p_nonwhite, na.rm = TRUE), 
            median_income = mean(median_income, na.rm = TRUE),
            total_population = mean(total_population, na.rm = TRUE))

call_geo = 
  calldata_merged_use %>%
  group_by(beat) %>%
  summarise(p_nonwhite = mean(p_nonwhite, na.rm = TRUE), 
            median_income = mean(median_income, na.rm = TRUE),
            total_population = mean(total_population, na.rm = TRUE))

## Separate by terciles
terry_geo$nw_high3 = 
  ifelse(terry_geo$p_nonwhite >= quantile(terry_geo$p_nonwhite, .67, na.rm = TRUE), 1, 0)
terry_geo$nw_low3 = 
  ifelse(terry_geo$p_nonwhite <= quantile(terry_geo$p_nonwhite, .33, na.rm = TRUE), 1, 0)

terry_geo$inc_high3 = 
  ifelse(terry_geo$median_income >= quantile(terry_geo$median_income, .67, na.rm = TRUE), 1, 0)
terry_geo$inc_low3 = 
  ifelse(terry_geo$median_income <= quantile(terry_geo$median_income, .33, na.rm = TRUE), 1, 0)

call_geo$nw_high3 = 
  ifelse(call_geo$p_nonwhite >= quantile(call_geo$p_nonwhite, .67, na.rm = TRUE), 1, 0)
call_geo$nw_low3 = 
  ifelse(call_geo$p_nonwhite <= quantile(call_geo$p_nonwhite, .33, na.rm = TRUE), 1, 0)

call_geo$inc_high3 = 
  ifelse(call_geo$median_income >= quantile(call_geo$median_income, .67, na.rm = TRUE), 1, 0)
call_geo$inc_low3 = 
  ifelse(call_geo$median_income <= quantile(call_geo$median_income, .33, na.rm = TRUE), 1, 0)

# Separate by median

terry_geo$nw_high2 = 
  ifelse(terry_geo$p_nonwhite >= quantile(terry_geo$p_nonwhite, .50, na.rm = TRUE), 1, 0)
terry_geo$nw_low2 = 
  ifelse(terry_geo$p_nonwhite <= quantile(terry_geo$p_nonwhite, .50, na.rm = TRUE), 1, 0)

terry_geo$inc_high2 = 
  ifelse(terry_geo$median_income >= quantile(terry_geo$median_income, .50, na.rm = TRUE), 1, 0)
terry_geo$inc_low2 = 
  ifelse(terry_geo$median_income <= quantile(terry_geo$median_income, .50, na.rm = TRUE), 1, 0)

call_geo$nw_high2 = 
  ifelse(call_geo$p_nonwhite >= quantile(call_geo$p_nonwhite, .50, na.rm = TRUE), 1, 0)
call_geo$nw_low2 = 
  ifelse(call_geo$p_nonwhite <= quantile(call_geo$p_nonwhite, .50, na.rm = TRUE), 1, 0)

call_geo$inc_high2 = 
  ifelse(call_geo$median_income >= quantile(call_geo$median_income, .50, na.rm = TRUE), 1, 0)
call_geo$inc_low2 = 
  ifelse(call_geo$median_income <= quantile(call_geo$median_income, .50, na.rm = TRUE), 1, 0)

# Separate geographic variables and merge with main data set

terry_geo = 
  terry_geo %>% 
  dplyr::select(nw_high3, nw_low3, inc_high3, inc_low3, nw_high2, nw_low2, inc_high2, inc_low2, beat)

call_geo = 
  call_geo %>% 
  dplyr::select(nw_high3, nw_low3, inc_high3, inc_low3, nw_high2, nw_low2, inc_high2, inc_low2, beat)

terry_geo_ts = terry_merged_use

call_geo_ts = calldata_merged_use

terry_geo_ts = merge(terry_geo_ts, terry_geo, by = "beat")
call_geo_ts = merge(call_geo_ts, call_geo, by = "beat")

## Separate terciles to calculate total pop within these

terry_geo3 = 
  terry_geo_ts %>%
  mutate(inc_3 = ifelse(inc_low3 == 1, "low", "high")) %>%
  mutate(nw_3 = ifelse(nw_low3 == 1, "low", "high"))

terry_geo_inc3 = 
  terry_geo3 %>%
  group_by(beat) %>% 
  summarise(total_population = mean(total_population))

terry_geo_inc3 = 
  terry_geo3 %>%
  group_by(inc_3) %>% 
  summarise(total_population_inc3 = mean(total_population, na.rm = TRUE))

terry_geo_nw3 = 
  terry_geo3 %>%
  group_by(nw_3) %>%
  summarise(total_population_nw3 = mean(total_population, na.rm = TRUE))

terry_geo3_ts = terry_geo3
terry_geo3_ts = merge(terry_geo3_ts, terry_geo_inc3, by = "inc_3")
terry_geo3_ts = merge(terry_geo3_ts, terry_geo_nw3, by = "nw_3")

call_geo3 = 
  call_geo_ts %>%
  mutate(inc_3 = ifelse(inc_low3 == 1, "low", "high")) %>%
  mutate(nw_3 = ifelse(nw_low3 == 1, "low", "high"))

call_geo_inc3 = 
  call_geo3 %>%
  group_by(inc_3) %>%
  summarise(total_population_inc3 = mean(total_population, na.rm = TRUE))

call_geo_nw3 = 
  call_geo3 %>%
  group_by(nw_3) %>%
  summarise(total_population_nw3 = mean(total_population, na.rm = TRUE))

call_geo3_ts = call_geo3
call_geo3_ts = merge(call_geo3_ts, call_geo_inc3, by = "inc_3")
call_geo3_ts = merge(call_geo3_ts, call_geo_nw3, by = "nw_3")

## Separate medians to calculate total pop within these

terry_geo2 = 
  terry_geo_ts %>%
  mutate(inc_2 = ifelse(inc_low2 == 1, "low", "high")) %>%
  mutate(nw_2 = ifelse(nw_low2 == 1, "low", "high"))

terry_geo_inc2 = 
  terry_geo2 %>%
  group_by(beat) %>% 
  summarise(total_population = mean(total_population))

terry_geo_inc2 = 
  terry_geo2 %>%
  group_by(inc_2) %>% 
  summarise(total_population_inc2 = mean(total_population, na.rm = TRUE))

terry_geo_nw2 = 
  terry_geo2 %>%
  group_by(nw_2) %>%
  summarise(total_population_nw2 = mean(total_population, na.rm = TRUE))

terry_geo2_ts = terry_geo2
terry_geo2_ts = merge(terry_geo2_ts, terry_geo_inc2, by = "inc_2")
terry_geo2_ts = merge(terry_geo2_ts, terry_geo_nw2, by = "nw_2")

call_geo2 = 
  call_geo_ts %>%
  mutate(inc_2 = ifelse(inc_low2 == 1, "low", "high")) %>%
  mutate(nw_2 = ifelse(nw_low2 == 1, "low", "high"))

call_geo_inc2 = 
  call_geo2 %>%
  group_by(inc_2) %>%
  summarise(total_population_inc2 = mean(total_population, na.rm = TRUE))

call_geo_nw2 = 
  call_geo2 %>%
  group_by(nw_2) %>%
  summarise(total_population_nw2 = mean(total_population, na.rm = TRUE))

call_geo2_ts = call_geo2
call_geo2_ts = merge(call_geo2_ts, call_geo_inc2, by = "inc_2")
call_geo2_ts = merge(call_geo2_ts, call_geo_nw2, by = "nw_2")

## Plot terciles

depol_geot_nwlow3 = terry_geo3_ts %>% 
  filter(nw_low3 == 1) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(total_population_nw3/1000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .5) + 
  labs(x = "Date", y = "Count",
       title = "A. Terry Stops")

depol_geot_nwhigh3 = terry_geo3_ts %>%
  filter(nw_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(total_population_nw3/1000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .5) + 
  labs(x = "Date", y = "Count",
       title = "A. Terry Stops")

depol_geoc_nwlow3 = call_geo3_ts %>% 
  filter(nw_low3 == 1) %>%
  mutate(count = 1) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(total_population_nw3/1000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             ize = .5) + 
  labs(x = "Date", y = "Count",
       title = "B. Officer-Initiated 911 Calls")

depol_geoc_nwhigh3 = call_geo3_ts %>% 
  filter(nw_high3 == 1) %>%
  mutate(count = 1) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(total_population_nw3/1000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             ize = .5) + 
  labs(x = "Date", y = "Count",
       title = "B. Officer-Initiated 911 Calls")

depol_geot_inclow3 = terry_geo3_ts %>% 
  filter(inc_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(total_population_inc3/1000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .5) + 
  labs(x = "Date", y = "Count",
       title = "A. Terry Stops")

depol_geot_inchigh3 = terry_geo3_ts %>% 
  filter(inc_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(total_population_inc3/1000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .5) + 
  labs(x = "Date", y = "Count",
       title = "A. Terry Stops")

depol_geoc_inclow3 = call_geo3_ts %>% 
  filter(inc_low3 == 1) %>%
  mutate(count = 1) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(total_population_inc3/1000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             ize = .5) + 
  labs(x = "Date", y = "Count",
       title = "B. Officer-Initiated 911 Calls")

depol_geoc_inchigh3 = call_geo3_ts %>% 
  filter(inc_high3 == 1) %>%
  mutate(count = 1) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(total_population_inc3/1000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             ize = .5) + 
  labs(x = "Date", y = "Count",
       title = "B. Officer-Initiated 911 Calls")

depol_nwlowGrob = arrangeGrob(depol_geot_nwlow3, depol_geoc_nwlow3, ncol = 2, top=textGrob("Percentage Nonwhite, Lowest Tercile"))
depol_nwhighGrob = arrangeGrob(depol_geot_nwhigh3, depol_geoc_nwhigh3, ncol = 2, top=textGrob("Percentage Nonwhite, Highest Tercile"))
depol_inclowGrob = arrangeGrob(depol_geot_inclow3, depol_geoc_inclow3, ncol = 2, top=textGrob("Median Income, Lowest Tercile"))
depol_inchighGrob = arrangeGrob(depol_geot_inchigh3, depol_geoc_inchigh3, ncol = 2, top=textGrob("Median Income, Highest Tercile"))

ggsave(plot = depol_nwlowGrob, filename = "depol_low3nonwhite.png", width = 8, height = 4)
ggsave(plot = depol_nwhighGrob, filename = "depol_high3nonwhite.png", width = 8, height = 4)
ggsave(plot = depol_inclowGrob, filename = "depol_low3income.png", width = 8, height = 4)
ggsave(plot = depol_inchighGrob, filename = "depol_high3income.png", width = 8, height = 4)

# Plot split by median

depol_geot_nwlow2 = terry_geo2_ts %>% 
  filter(nw_low2 == 1) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(total_population_nw2/1000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
             #ylim(0,0.004) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .5) + 
  labs(x = "Date", y = "Count",
       title = "A. Terry Stops")

depol_geot_nwhigh2 = terry_geo2_ts %>% 
  filter(nw_high2 == 1) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(total_population_nw2/1000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
             #ylim(0,0.004) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .5) + 
  labs(x = "Date", y = "Count",
       title = "A. Terry Stops")

depol_geot_inclow2 = terry_geo2_ts %>% 
  filter(inc_low2 == 1) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(total_population_inc2/1000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
             #ylim(0,0.004) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .5) + 
  labs(x = "Date", y = "Count",
       title = "A. Terry Stops")

depol_geot_inchigh2 = terry_geo2_ts %>% 
  filter(inc_high2 == 1) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(total_population_inc2/1000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
             #ylim(0,0.004) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .5) + 
  labs(x = "Date", y = "Count",
       title = "A. Terry Stops")

depol_geoc_nwlow2 = call_geo2_ts %>% 
  filter(nw_low2 == 1) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(total_population_nw2/1000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
             #ylim(0,0.004) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .5) + 
  labs(x = "Date", y = "Count",
       title = "B. Officer-Initiated 911 Calls")

depol_geoc_nwhigh2 = call_geo2_ts %>% 
  filter(nw_high2 == 1) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(total_population_nw2/1000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
             #ylim(0,0.004) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .5) + 
  labs(x = "Date", y = "Count",
       title = "B. Officer-Initiated 911 Calls")

depol_geoc_inclow2 = call_geo2_ts %>% 
  filter(inc_low2 == 1) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(total_population_inc2/1000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
             #ylim(0,0.004) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .5) + 
  labs(x = "Date", y = "Count",
       title = "B. Officer-Initiated 911 Calls")

depol_geoc_inchigh2 = call_geo2_ts %>% 
  filter(inc_high2 == 1) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(total_population_inc2/1000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
             #ylim(0,0.004) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .5) + 
  labs(x = "Date", y = "Count",
       title = "B. Officer-Initiated 911 Calls")

## Arrange 

depol_nwlow2Grob = arrangeGrob(depol_geot_nwlow2, depol_geoc_nwlow2, ncol = 2, top=textGrob("Percentage Nonwhite, Below Median"))
depol_nwhigh2Grob = arrangeGrob(depol_geot_nwhigh2, depol_geoc_nwhigh2, ncol = 2, top=textGrob("Percentage Nonwhite, Above Median"))
depol_inclow2Grob = arrangeGrob(depol_geot_inclow2, depol_geoc_inclow2, ncol = 2, top=textGrob("Median Income, Below Median"))
depol_inchigh2Grob = arrangeGrob(depol_geot_inchigh2, depol_geoc_inchigh2, ncol = 2, top=textGrob("Median Income, Above Median"))

ggsave(plot = depol_nwlow2Grob, filename = "depol_low2nonwhite.png", width = 8, height = 4)
ggsave(plot = depol_nwhigh2Grob, filename = "depol_high2nonwhite.png", width = 8, height = 4)
ggsave(plot = depol_inclow2Grob, filename = "depol_low2income.png", width = 8, height = 4)
ggsave(plot = depol_inchigh2Grob, filename = "depol_high2income.png", width = 8, height = 4)

#### Analysis: RDs 25, 50, 100 bw
terry_geo2_ts_25 = terry_geo2_ts %>% 
  filter(date >= as.Date("2020-05-04") & date <= as.Date("2020-06-22"))

terry_geo2_ts_50 = terry_geo2_ts %>% 
  filter(date >= as.Date("2020-04-08") & date <= as.Date("2020-07-17"))

terry_geo2_ts_100 = terry_geo2_ts %>% 
  filter(date >= as.Date("2020-02-17") & date <= as.Date("2020-09-05"))

terry_geo2_nwlow_25 = terry_geo2_ts %>% 
  filter(date >= as.Date("2020-05-04") & date <= as.Date("2020-06-22")) %>% 
  filter(nw_low2 == 1)

terry_geo2_nwlow_50 = terry_geo2_ts %>% 
  filter(date >= as.Date("2020-04-08") & date <= as.Date("2020-07-17")) %>% 
  filter(nw_low2 == 1)

terry_geo2_nwlow_100 = terry_geo2_ts %>% 
  filter(date >= as.Date("2020-02-17") & date <= as.Date("2020-09-05")) %>% 
  filter(nw_low2 == 1)

terry_geo2_nwhigh_25 = terry_geo2_ts %>% 
  filter(date >= as.Date("2020-05-04") & date <= as.Date("2020-06-22")) %>% 
  filter(nw_high2 == 1)

terry_geo2_nwhigh_50 = terry_geo2_ts %>% 
  filter(date >= as.Date("2020-04-08") & date <= as.Date("2020-07-17")) %>% 
  filter(nw_high2 == 1)

terry_geo2_nwhigh_100 = terry_geo2_ts %>% 
  filter(date >= as.Date("2020-02-17") & date <= as.Date("2020-09-05")) %>% 
  filter(nw_high2 == 1)

terry_geo2_inclow_25 = terry_geo2_ts %>% 
  filter(date >= as.Date("2020-05-04") & date <= as.Date("2020-06-22")) %>% 
  filter(inc_low2 == 1)

terry_geo2_inclow_50 = terry_geo2_ts %>% 
  filter(date >= as.Date("2020-04-08") & date <= as.Date("2020-07-17")) %>% 
  filter(inc_low2 == 1)

terry_geo2_inclow_100 = terry_geo2_ts %>% 
  filter(date >= as.Date("2020-02-17") & date <= as.Date("2020-09-05")) %>% 
  filter(inc_low2 == 1)

terry_geo2_inchigh_25 = terry_geo2_ts %>% 
  filter(date >= as.Date("2020-05-04") & date <= as.Date("2020-06-22")) %>% 
  filter(inc_high2 == 1)

terry_geo2_inchigh_50 = terry_geo2_ts %>% 
  filter(date >= as.Date("2020-04-08") & date <= as.Date("2020-07-17")) %>% 
  filter(inc_high2 == 1)

terry_geo2_inchigh_100 = terry_geo2_ts %>% 
  filter(date >= as.Date("2020-02-17") & date <= as.Date("2020-09-05")) %>% 
  filter(inc_high2 == 1)

call_geo2_ts_25 = call_geo2_ts %>% 
  filter(date >= as.Date("2020-05-04") & date <= as.Date("2020-06-22"))

call_geo2_ts_50 = call_geo2_ts %>% 
  filter(date >= as.Date("2020-04-08") & date <= as.Date("2020-07-17"))

call_geo2_ts_100 = call_geo2_ts %>% 
  filter(date >= as.Date("2020-02-17") & date <= as.Date("2020-09-05"))

call_geo2_nwlow_25 = call_geo2_ts %>% 
  filter(date >= as.Date("2020-05-04") & date <= as.Date("2020-06-22")) %>% 
  filter(nw_low2 == 1)

call_geo2_nwlow_50 = call_geo2_ts %>% 
  filter(date >= as.Date("2020-04-08") & date <= as.Date("2020-07-17")) %>% 
  filter(nw_low2 == 1)

call_geo2_nwlow_100 = call_geo2_ts %>% 
  filter(date >= as.Date("2020-02-17") & date <= as.Date("2020-09-05")) %>% 
  filter(nw_low2 == 1)

call_geo2_nwhigh_25 = call_geo2_ts %>% 
  filter(date >= as.Date("2020-05-04") & date <= as.Date("2020-06-22")) %>% 
  filter(nw_high2 == 1)

call_geo2_nwhigh_50 = call_geo2_ts %>% 
  filter(date >= as.Date("2020-04-08") & date <= as.Date("2020-07-17")) %>% 
  filter(nw_high2 == 1)

call_geo2_nwhigh_100 = call_geo2_ts %>% 
  filter(date >= as.Date("2020-02-17") & date <= as.Date("2020-09-05")) %>% 
  filter(nw_high2 == 1)

call_geo2_inclow_25 = call_geo2_ts %>% 
  filter(date >= as.Date("2020-05-04") & date <= as.Date("2020-06-22")) %>% 
  filter(inc_low2 == 1)

call_geo2_inclow_50 = call_geo2_ts %>% 
  filter(date >= as.Date("2020-04-08") & date <= as.Date("2020-07-17")) %>% 
  filter(inc_low2 == 1)

call_geo2_inclow_100 = call_geo2_ts %>% 
  filter(date >= as.Date("2020-02-17") & date <= as.Date("2020-09-05")) %>% 
  filter(inc_low2 == 1)

call_geo2_inchigh_25 = call_geo2_ts %>% 
  filter(date >= as.Date("2020-05-04") & date <= as.Date("2020-06-22")) %>% 
  filter(inc_high2 == 1)

call_geo2_inchigh_50 = call_geo2_ts %>% 
  filter(date >= as.Date("2020-04-08") & date <= as.Date("2020-07-17")) %>% 
  filter(inc_high2 == 1)

call_geo2_inchigh_100 = call_geo2_ts %>% 
  filter(date >= as.Date("2020-02-17") & date <= as.Date("2020-09-05")) %>% 
  filter(inc_high2 == 1)

## Separate by outcomes
terry_geo_inc2_low_ts = terry_geo2_ts %>% 
  filter(inc_low2 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(total_population_inc2/1000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

terry_geo_inc2_high_ts = terry_geo2_ts %>% 
  filter(inc_high2 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(total_population_inc2/1000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

terry_geo_nw2_low_ts = terry_geo2_ts %>% 
  filter(nw_low2 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(total_population_nw2/1000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

terry_geo_nw2_high_ts = terry_geo2_ts %>% 
  filter(nw_high2 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(total_population_nw2/1000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

call_geo_inc2_low_ts = call_geo2_ts %>% 
  filter(inc_low2 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(total_population_inc2/1000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

call_geo_inc2_high_ts = call_geo2_ts %>% 
  filter(inc_high2 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(total_population_inc2/1000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

call_geo_nw2_low_ts = call_geo2_ts %>% 
  filter(nw_low2 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(total_population_nw2/1000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

call_geo_nw2_high_ts = call_geo2_ts %>% 
  filter(nw_high2 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(total_population_nw2/1000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

terry_geo_inc2_low_ts$blmrv = 
  terry_geo_inc2_low_ts$date - min(terry_geo_inc2_low_ts$date[terry_geo_inc2_low_ts$blm == 1])
terry_geo_inc2_high_ts$blmrv = 
  terry_geo_inc2_high_ts$date - min(terry_geo_inc2_high_ts$date[terry_geo_inc2_high_ts$blm == 1])
terry_geo_nw2_low_ts$blmrv = 
  terry_geo_nw2_low_ts$date - min(terry_geo_nw2_low_ts$date[terry_geo_nw2_low_ts$blm == 1])
terry_geo_nw2_high_ts$blmrv = 
  terry_geo_nw2_high_ts$date - min(terry_geo_nw2_high_ts$date[terry_geo_nw2_high_ts$blm == 1])

call_geo_inc2_low_ts$blmrv = 
  call_geo_inc2_low_ts$date - min(call_geo_inc2_low_ts$date[call_geo_inc2_low_ts$blm == 1])
call_geo_inc2_high_ts$blmrv = 
  call_geo_inc2_high_ts$date - min(call_geo_inc2_high_ts$date[call_geo_inc2_high_ts$blm == 1])
call_geo_nw2_low_ts$blmrv = 
  call_geo_nw2_low_ts$date - min(call_geo_nw2_low_ts$date[call_geo_nw2_low_ts$blm == 1])
call_geo_nw2_high_ts$blmrv = 
  call_geo_nw2_high_ts$date - min(call_geo_nw2_high_ts$date[call_geo_nw2_high_ts$blm == 1])

datasets_depol_geo = list(terry_geo_inc2_low_ts %>% as.data.frame %>% mutate(count_terry_inc_low = count) %>% filter(date >= as.Date("2016-01-01")),
                          terry_geo_inc2_high_ts %>% as.data.frame %>% mutate(count_terry_inc_high = count) %>% filter(date >= as.Date("2016-01-01")),
                          terry_geo_nw2_low_ts %>% as.data.frame %>% mutate(count_terry_nw_low = count) %>% filter(date >= as.Date("2016-01-01")),
                          terry_geo_nw2_high_ts %>% as.data.frame %>% mutate(count_terry_nw_high = count) %>% filter(date >= as.Date("2016-01-01")),
                          call_geo_inc2_low_ts %>% as.data.frame %>% mutate(count_calls_inc_low = count) %>% filter(date >= as.Date("2016-01-01")),
                          call_geo_inc2_high_ts %>% as.data.frame %>% mutate(count_calls_inc_high = count) %>% filter(date >= as.Date("2016-01-01")),
                          call_geo_nw2_low_ts %>% as.data.frame %>% mutate(count_calls_nw_low = count) %>% filter(date >= as.Date("2016-01-01")),
                          call_geo_nw2_high_ts %>% as.data.frame %>% mutate(count_calls_nw_high = count) %>% filter(date >= as.Date("2016-01-01")))

# outcomes

outcomes = c("count_terry_inc_low", "count_terry_inc_high", 
             "count_terry_nw_low", "count_terry_nw_high",
             "count_calls_inc_low", "count_calls_inc_high",
             "count_calls_nw_low", "count_calls_nw_high")

# polynomials

polys = c(0, 1, 2)

# kernels 

kerns = c("tri", "uni", "epa")

# loop for analysis 

dataset_list = as.list(rep(NA, length(datasets_depol_geo)))

for (i in 1:length(datasets_depol_geo)) {
  
  poly_list = as.list(rep(NA, length(polys)))
  
  for (j in 1:length(polys)) {
    
    print(paste0("J Iteration ", j))
    
    kern_list = as.list(rep(NA, length(kerns)))
    
    for (l in 1:length(kerns)) {
      
      print(paste0("L Iteration ", l))
      
      mat_bw = matrix(NA, 100, 7)
      
      for(z in c(25, 50, 100)) {
        
        print(paste0("Z Iteration ", z))
        
        outrd = 
          rdrobust(x = as.numeric(datasets_depol_geo[[i]][, "blmrv"]),
                   y = datasets_depol_geo[[i]][, outcomes[i]],
                   p = polys[j],
                   h = z,
                   kernel = kerns[l])
        
        mat_bw[z, 1] = outrd$coef[3]
        mat_bw[z, 2] = outrd$se[3]
        mat_bw[z, 3] = outrd$pv[3]
        mat_bw[z, 4] = polys[j]
        mat_bw[z, 5] = kerns[l]
        mat_bw[z, 6] = outcomes[i]
        mat_bw[z, 7] = z
        
      }
      
      kern_list[[l]] = mat_bw %>%
        as.data.frame() %>%
        `colnames<-` (c("Coefficient", "SE", "Pvalue", "Polynomial", "Kernel", "Outcome", "BW"))%>% 
        mutate(Coefficient = as.numeric(as.character(Coefficient)),
               SE = as.numeric(as.character(SE)),
               Pvalue = as.numeric(as.character(Pvalue)),
               Polynomial = as.numeric(as.character(Polynomial)),
               Kernel = as.character(Kernel),
               Outcome = as.character(Outcome),
               BW = as.character(BW))
    }
    
    poly_list[[j]] = kern_list
    
  }
  
  dataset_list[[i]] = poly_list
  
}

depol_rdit_geo = dataset_list %>% 
  unlist(recursive = FALSE) %>%
  unlist(recursive = FALSE) %>%
  do.call(rbind.data.frame, .) %>% 
  mutate(Outcome = factor(Outcome,
                          levels = c("count_terry_inc_low", "count_terry_inc_high", 
                                     "count_terry_nw_low", "count_terry_nw_high",
                                     "count_calls_inc_low", "count_calls_inc_high",
                                     "count_calls_nw_low", "count_calls_nw_high"),
                          labels = c("A. Terry Stops, Low Income",
                                     "B. Terry Stops, High Income",
                                     "C. Terry Stops, Low Percentage Nonwhite",
                                     "D. Terry Stops, High Percentage Nonwhite",
                                     "E. Officer Initiated 911 Calls, Low Income",
                                     "F. Officer Initiated 911 Calls, High Income",
                                     "G. Officer Initiated 911 Calls, Low Percentage Nonwhite",
                                     "H. Officer Initiated 911 Calls, High Percentage Nonwhite"
                                     )),
         Polynomial = factor(Polynomial, levels = c(0, 1, 2),
                             labels = c("DIM", "Linear", "Quadratic")),
         Kernel = factor(Kernel, levels = c("tri", "uni", "epa"),
                         labels = c("Triangular",
                                    "Uniform",
                                    "Epanechnikov"))) %>% 
  filter(!is.na(BW))

depol_rdit_geo = depol_rdit_geo %>% 
  mutate(kern_poly = factor(paste0(Kernel, ", ", Polynomial), levels = c("Triangular, DIM", 
                                                                         "Triangular, Linear",
                                                                         "Triangular, Quadratic",
                                                                         "Uniform, DIM",
                                                                         "Uniform, Linear",
                                                                         "Uniform, Quadratic",
                                                                         "Epanechnikov, DIM",
                                                                         "Epanechnikov, Linear",
                                                                         "Epanechnikov, Quadratic")))

depol_rdit_geo$BW = as.numeric(depol_rdit_geo$BW)

#### Coefficient tests of difference ####
## Create function for s
s <- function(n1, n2, p, s1, s2){sqrt(((n1 - p)*((s1)^2) + (n2- p)*((s2)^2)/(n1 + n2 - 2*p)))}

se <- function(s, SEb1, SEb2, s1, s2){s*sqrt(((SEb1/s1)^2) + (SEb2/s2)^2)}

z <- function(b1, b2, se){(b1-b2)/se}

b_est <- function(b1, b2){(b1-b2)}

s_terry_25_inc <- s(n1 = 467, n2 = 359, p = 2, s1 = 0.4583717, s2 = 0.2447864)
se_t_25_inc <- se(s = s_terry_25_inc, SEb1 = 0.4583717, SEb2 = 0.2447864, s1 = 0.4583717, s2 = 0.2447864)
z_t_25_inc <- z(b1 = -1.21717345, b2 = -0.33903884, se = se_t_25_inc)
b_t_25_inc <- b_est(b1 = -1.21717345, b2 = -0.33903884)
s_terry_25_nw <- s(n1 = 375, n2 = 451, p = 2, s1 = 0.2166516, s2 = 0.5313162)
se_t_25_nw <- se(s = s_terry_25_nw, SEb1 = 0.2166516, SEb2 = 0.5313162, s1 = 0.2166516, s2 = 0.5313162)
z_t_25_nw <- z(b1 = -0.18904011, b2 = -1.64855498, se = se_t_25_nw)
b_t_25_nw <- b_est(b1 = -0.18904011, b2 = -1.64855498)

s_terry_50_inc <- s(n1 = 1022, n2 = 774, p = 2, s1 = 0.3187507, s2 = 0.1664611)
se_t_50_inc <- se(s = s_terry_50_inc, SEb1 = 0.3187507, SEb2 = 0.1664611, s1 = 0.3187507, s2 = 0.1664611)
z_t_50_inc <- z(b1 = -1.30696829, b2 = -0.33307122, se = se_t_50_inc)
b_t_50_inc <- b_est(b1 = -1.30696829, b2 = -0.33307122)
s_terry_50_nw <- s(n1 = 783, n2 = 1013, p = 2, s1 = 0.1532174, s2 = 0.3506309)
se_t_50_nw <- se(s = s_terry_50_nw, SEb1 = 0.1532174, SEb2 = 0.3506309, s1 = 0.1532174, s2 = 0.3506309)
z_t_50_nw <- z(b1 = -0.34365118, b2 = -1.41273639, se = se_t_50_nw)
b_t_50_nw <- b_est(b1 = -0.34365118, b2 = -1.41273639)

s_terry_100_inc <- s(n1 = 1881, n2 = 1489, p = 2, s1 = 0.2387584, s2 = 0.1152557)
se_t_100_inc <- se(s = s_terry_100_inc, SEb1 = 0.2387584, SEb2 = 0.1152557, s1 = 0.2387584, s2 = 0.1152557)
z_t_100_inc <- z(b1 = -1.80616589, b2 = -0.59693420, se = se_t_100_inc)
b_t_100_inc <- b_est(b1 = -1.80616589, b2 = -0.59693420)
s_terry_100_nw <- s(n1 = 1439, n2 = 1931, p = 2, s1 = 0.1170986, s2 = 0.2493603)
se_t_100_nw <- se(s = s_terry_100_nw, SEb1 = 0.1170986, SEb2 = 0.2493603, s1 = 0.1170986, s2 = 0.2493603)
z_t_100_nw <- z(b1 = -0.64595796, b2 = -1.81168528, se = se_t_100_nw)
b_t_100_nw <- b_est(b1 = -0.64595796, b2 = -1.81168528)

s_call_25_inc <- s(n1 = 9078, n2 = 5948, p = 2, s1 = 6.2913145, s2 = 3.3857775)
se_c_25_inc <- se(s = s_call_25_inc, SEb1 = 6.2913145, SEb2 = 3.3857775, s1 = 6.2913145, s2 = 3.3857775)
z_c_25_inc <- z(b1 = -26.27996963, b2 = -10.17270712, se = se_c_25_inc)
b_c_25_inc <- b_est(b1 = -26.27996963, b2 = -10.17270712)
s_call_25_nw <- s(n1 = 6043, n2 = 8983, p = 2, s1 = 4.2101540, s2 = 4.5415662)
se_c_25_nw <- se(s = s_call_25_nw, SEb1 = 4.2101540, SEb2 = 4.5415662, s1 = 4.2101540, s2 = 4.5415662)
z_c_25_nw <- z(b1 = -10.50082801, b2 = -24.44887013, se = se_c_25_nw)
b_c_25_nw <- b_est(b1 = -10.50082801, b2 = -24.44887013)

s_call_50_inc <- s(n1 = 19040, n2 = 12767, p = 2, s1 = 3.8690116, s2 = 1.8795014)
se_c_50_inc <- se(s = s_call_50_inc, SEb1 = 3.8690116, SEb2 = 1.8795014, s1 = 3.8690116, s2 = 1.8795014)
z_c_50_inc <- z(b1 = -32.58404515, b2 = -12.95343923, se = se_c_50_inc)
b_c_50_inc <- b_est(b1 = -32.58404515, b2 = -12.95343923)
s_call_50_nw <- s(n1 = 12923, n2 = 18884, p = 2, s1 = 2.3495744, s2 = 2.8426370)
se_c_50_nw <- se(s = s_call_50_nw, SEb1 = 2.3495744, SEb2 = 2.8426370, s1 = 2.3495744, s2 = 2.8426370)
z_c_50_nw <- z(b1 = -13.25036853, b2 = -30.52253813, se = se_c_50_nw)
b_c_50_nw <- b_est(b1 = -13.25036853, b2 = -30.52253813)

s_call_100_inc <- s(n1 = 38941, n2 = 25054, p = 2, s1 = 2.6143121, s2 = 1.0786114)
se_c_100_inc <- se(s = s_call_100_inc, SEb1 = 2.6143121, SEb2 = 1.0786114, s1 = 2.6143121, s2 = 1.0786114)
z_c_100_inc <- z(b1 = -37.12474311, b2 = -15.09381730, se = se_c_100_inc)
b_c_100_inc <- b_est(b1 = -37.12474311, b2 = -15.09381730)
s_call_100_nw <- s(n1 = 24507, n2 = 39488, p = 2, s1 = 1.3538813, s2 = 2.0196423)
se_c_100_nw <- se(s = s_call_100_nw, SEb1 = 1.3538813, SEb2 = 2.0196423, s1 = 1.3538813, s2 = 2.0196423)
z_c_100_nw <- z(b1 = -15.14702069, b2 = -35.27757777, se = se_c_100_nw)
b_c_100_nw <- b_est(b1 = -15.14702069, b2 = -35.27757777)


#### produce bivariate plots for paper ####

# calls
biv_call = call %>% 
  filter(initiate_off == 1) %>%
  filter(date >= as.Date("2016-01-01")) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .5) + 
  labs(x = "Date", y = "Count",
       title = "A. Officer-Initiated 911 Calls")

# stops
biv_terry = terry %>% 
  mutate(count = 1) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .5) + 
  labs(x = "Date", y = "Count",
       title = "B. Terry Stops")

biv_depolGrob = arrangeGrob(biv_call, biv_terry, ncol = 2)
ggsave(plot = biv_depolGrob, filename = "biv_depol.png", width = 8, height = 3)

# arrest rate
biv_arrest = terry %>% 
  filter(date >= as.Date("2019-01-01")) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            arrest = mean(arrest, na.rm = TRUE)) %>% 
  ggplot() + 
  geom_point(aes(x = date, y = arrest),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = arrest, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .5) + 
  labs(x = "Date", y = "Arrest Rate",
       title = "A. Arrest Rate")

# hit rate
biv_hit = terry %>% 
  filter(date >= as.Date("2019-01-01")) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            hit = mean(hit, na.rm = TRUE)) %>% 
  ggplot() + 
  geom_point(aes(x = date, y = hit),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = hit, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .5) + 
  labs(x = "Date", y = "Hit Rate",
       title = "B. Hit Rate")

# bw disp
biv_disp = terry_arr_ts %>% 
  filter(date >= as.Date("2019-01-01")) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            rr_bw = mean(rr_bw, na.rm = TRUE)) %>% 
  ggplot() + 
  geom_point(aes(x = date, y = rr_bw),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = rr_bw, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .5) + 
  labs(x = "Date", y = "Rate Ratio",
       title = "C. Black/White Disparity")


# resp time
biv_resp = call %>% 
  filter(date >= as.Date("2019-01-01")) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            time = mean(response_time, na.rm = TRUE)) %>% 
  ggplot() + 
  geom_point(aes(x = date, y = time),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = time, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .5) + 
  labs(x = "Date", y = "Seconds",
       title = "D. Response Time")

biv_effGrob = arrangeGrob(biv_arrest, biv_hit, biv_disp, biv_resp, ncol = 2)
ggsave(plot = biv_effGrob, filename = "biv_eff.png", width = 8, height = 6)

# all crime
biv_crm = crime %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .5) + 
  labs(x = "Date", y = "Count",
       title = "A. All Crimes")

# against person crime
biv_pers_crm = crime %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  filter(cat_person == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .5) + 
  labs(x = "Date", y = "Count",
       title = "B. Against Person Crimes")

# against property crime
biv_prop_crm = crime %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  filter(cat_property == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .5) + 
  labs(x = "Date", y = "Count",
       title = "C. Against Property Crimes")

# against society crime
biv_soc_crm = crime %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  filter(cat_society == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .5) + 
  labs(x = "Date", y = "Count",
       title = "D. Against Society Crimes")

biv_crmGrob = arrangeGrob(biv_crm, biv_pers_crm, biv_prop_crm, biv_soc_crm, ncol = 2)
ggsave(plot = biv_crmGrob, filename = "biv_crm.png", width = 8, height = 6)
